/*This file is part of the FEBio source code and is licensed under the MIT license
listed below.

See Copyright-FEBio.txt for details.

Copyright (c) 2020 University of Utah, The Trustees of Columbia University in 
the City of New York, and others.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.*/



#pragma once
// NOTE: This file is automatically included from tens4d.h
// Users should not include this file manually!

inline tens4d::tens4d(const double g)
{
    d[ 0] = d[ 1] = d[ 2] = d[ 3] = d[ 4] = d[ 5] = d[ 6] = d[ 7] = d[ 8] = g;
    d[ 9] = d[10] = d[11] = d[12] = d[13] = d[14] = d[15] = d[16] = d[17] = g;
    d[18] = d[19] = d[20] = d[21] = d[22] = d[23] = d[24] = d[25] = d[26] = g;
    d[27] = d[28] = d[29] = d[30] = d[31] = d[32] = d[33] = d[34] = d[35] = g;
    d[36] = d[37] = d[38] = d[39] = d[40] = d[41] = d[42] = d[43] = d[44] = g;
    d[45] = d[46] = d[47] = d[48] = d[49] = d[50] = d[51] = d[52] = d[53] = g;
    d[54] = d[55] = d[56] = d[57] = d[58] = d[59] = d[60] = d[61] = d[62] = g;
    d[63] = d[64] = d[65] = d[66] = d[67] = d[68] = d[69] = d[70] = d[71] = g;
    d[72] = d[73] = d[74] = d[75] = d[76] = d[77] = d[78] = d[79] = d[80] = g;
}

inline tens4d::tens4d(const tens4ds t)
{
    d[ 0]=t.d[ 0]; d[ 9]=t.d[ 1]; d[18]=t.d[ 3]; d[27]=t.d[ 6]; d[36]=t.d[10]; d[45]=t.d[15]; d[54]=t.d[ 6]; d[63]=t.d[10]; d[72]=t.d[15];
    d[ 1]=t.d[ 1]; d[10]=t.d[ 2]; d[19]=t.d[ 4]; d[28]=t.d[ 7]; d[37]=t.d[11]; d[46]=t.d[16]; d[55]=t.d[ 7]; d[64]=t.d[11]; d[73]=t.d[16];
    d[ 2]=t.d[ 3]; d[11]=t.d[ 4]; d[20]=t.d[ 5]; d[29]=t.d[ 8]; d[38]=t.d[12]; d[47]=t.d[17]; d[56]=t.d[ 8]; d[65]=t.d[12]; d[74]=t.d[17];
    d[ 3]=t.d[ 6]; d[12]=t.d[ 7]; d[21]=t.d[ 8]; d[30]=t.d[ 9]; d[39]=t.d[13]; d[48]=t.d[18]; d[57]=t.d[ 9]; d[66]=t.d[13]; d[75]=t.d[18];
    d[ 4]=t.d[10]; d[13]=t.d[11]; d[22]=t.d[12]; d[31]=t.d[13]; d[40]=t.d[14]; d[49]=t.d[19]; d[58]=t.d[13]; d[67]=t.d[14]; d[76]=t.d[19];
    d[ 5]=t.d[15]; d[14]=t.d[16]; d[23]=t.d[17]; d[32]=t.d[18]; d[41]=t.d[19]; d[50]=t.d[20]; d[59]=t.d[18]; d[68]=t.d[19]; d[77]=t.d[20];
    d[ 6]=t.d[ 6]; d[15]=t.d[ 7]; d[24]=t.d[ 8]; d[33]=t.d[ 9]; d[42]=t.d[13]; d[51]=t.d[18]; d[60]=t.d[ 9]; d[69]=t.d[13]; d[78]=t.d[18];
    d[ 7]=t.d[10]; d[16]=t.d[11]; d[25]=t.d[12]; d[34]=t.d[13]; d[43]=t.d[14]; d[52]=t.d[19]; d[61]=t.d[13]; d[70]=t.d[14]; d[79]=t.d[19];
    d[ 8]=t.d[15]; d[17]=t.d[16]; d[26]=t.d[17]; d[35]=t.d[18]; d[44]=t.d[19]; d[53]=t.d[20]; d[62]=t.d[18]; d[71]=t.d[19]; d[80]=t.d[20];
}

inline tens4d::tens4d(const tens4dmm t)
{
    d[ 0]=t.d[ 0]; d[ 9]=t.d[ 6]; d[18]=t.d[12]; d[27]=t.d[18]; d[36]=t.d[24]; d[45]=t.d[30]; d[54]=t.d[18]; d[63]=t.d[24]; d[72]=t.d[30];
    d[ 1]=t.d[ 1]; d[10]=t.d[ 7]; d[19]=t.d[13]; d[28]=t.d[19]; d[37]=t.d[25]; d[46]=t.d[31]; d[55]=t.d[19]; d[64]=t.d[25]; d[73]=t.d[31];
    d[ 2]=t.d[ 2]; d[11]=t.d[ 8]; d[20]=t.d[14]; d[29]=t.d[20]; d[38]=t.d[26]; d[47]=t.d[32]; d[56]=t.d[20]; d[65]=t.d[26]; d[74]=t.d[32];
    d[ 3]=t.d[ 3]; d[12]=t.d[ 9]; d[21]=t.d[15]; d[30]=t.d[21]; d[39]=t.d[27]; d[48]=t.d[33]; d[57]=t.d[21]; d[66]=t.d[27]; d[75]=t.d[33];
    d[ 4]=t.d[ 4]; d[13]=t.d[10]; d[22]=t.d[16]; d[31]=t.d[22]; d[40]=t.d[28]; d[49]=t.d[34]; d[58]=t.d[22]; d[67]=t.d[28]; d[76]=t.d[34];
    d[ 5]=t.d[ 5]; d[14]=t.d[11]; d[23]=t.d[17]; d[32]=t.d[23]; d[41]=t.d[29]; d[50]=t.d[35]; d[59]=t.d[23]; d[68]=t.d[29]; d[77]=t.d[35];
    d[ 6]=t.d[ 3]; d[15]=t.d[ 9]; d[24]=t.d[15]; d[33]=t.d[21]; d[42]=t.d[27]; d[51]=t.d[33]; d[60]=t.d[21]; d[69]=t.d[27]; d[78]=t.d[33];
    d[ 7]=t.d[ 4]; d[16]=t.d[10]; d[25]=t.d[16]; d[34]=t.d[22]; d[43]=t.d[28]; d[52]=t.d[34]; d[61]=t.d[22]; d[70]=t.d[28]; d[79]=t.d[34];
    d[ 8]=t.d[ 5]; d[17]=t.d[11]; d[26]=t.d[17]; d[35]=t.d[23]; d[44]=t.d[29]; d[53]=t.d[35]; d[62]=t.d[23]; d[71]=t.d[29]; d[80]=t.d[35];
}

inline tens4d::tens4d(double m[9][9])
{
    d[ 0]=m[0][0]; d[ 9]=m[0][1]; d[18]=m[0][2]; d[27]=m[0][3]; d[36]=m[0][4]; d[45]=m[0][5]; d[54]=m[0][6]; d[63]=m[0][7]; d[72]=m[0][8];
    d[ 1]=m[1][0]; d[10]=m[1][1]; d[19]=m[1][2]; d[28]=m[1][3]; d[37]=m[1][4]; d[46]=m[1][5]; d[55]=m[1][6]; d[64]=m[1][7]; d[73]=m[1][8];
    d[ 2]=m[2][0]; d[11]=m[2][1]; d[20]=m[2][2]; d[29]=m[2][3]; d[38]=m[2][4]; d[47]=m[2][5]; d[56]=m[2][6]; d[65]=m[2][7]; d[74]=m[2][8];
    d[ 3]=m[3][0]; d[12]=m[3][1]; d[21]=m[3][2]; d[30]=m[3][3]; d[39]=m[3][4]; d[48]=m[3][5]; d[57]=m[3][6]; d[66]=m[3][7]; d[75]=m[3][8];
    d[ 4]=m[4][0]; d[13]=m[4][1]; d[22]=m[4][2]; d[31]=m[4][3]; d[40]=m[4][4]; d[49]=m[4][5]; d[58]=m[4][6]; d[67]=m[4][7]; d[76]=m[4][8];
    d[ 5]=m[5][0]; d[14]=m[5][1]; d[23]=m[5][2]; d[32]=m[5][3]; d[41]=m[5][4]; d[50]=m[5][5]; d[59]=m[5][6]; d[68]=m[5][7]; d[77]=m[5][8];
    d[ 6]=m[6][0]; d[15]=m[6][1]; d[24]=m[6][2]; d[33]=m[6][3]; d[42]=m[6][4]; d[51]=m[6][5]; d[60]=m[6][6]; d[69]=m[6][7]; d[78]=m[6][8];
    d[ 7]=m[7][0]; d[16]=m[7][1]; d[25]=m[7][2]; d[34]=m[7][3]; d[43]=m[7][4]; d[52]=m[7][5]; d[61]=m[7][6]; d[70]=m[7][7]; d[79]=m[7][8];
    d[ 8]=m[8][0]; d[17]=m[8][1]; d[26]=m[8][2]; d[35]=m[8][3]; d[44]=m[8][4]; d[53]=m[8][5]; d[62]=m[8][6]; d[71]=m[8][7]; d[80]=m[8][8];
}

inline double& tens4d::operator () (int i, int j, int k, int l)
{
	const int m[3][3] = {{0,3,8},{6,1,4},{5,7,2}};
	tens4d& T = (*this);
	return T(m[i][j], m[k][l]);
}

inline double tens4d::operator () (int i, int j, int k, int l) const
{
	const int m[3][3] = {{0,3,8},{6,1,4},{5,7,2}};
	const tens4d& T = (*this);
	return T(m[i][j], m[k][l]);
}

inline double& tens4d::operator () (int i, int j)
{
	const int m[9] = {0, 9, 18, 27, 36, 45, 54, 63, 72};
	return d[m[j]+i];
}

inline double tens4d::operator () (int i, int j) const
{
	const int m[9] = {0, 9, 18, 27, 36, 45, 54, 63, 72};
	return d[m[j]+i];
}

//-----------------------------------------------------------------------------
// operator +
inline tens4d tens4d::operator + (const tens4d& t) const
{
    tens4d s;
//    for (int i=0; i<NNZ; i++)
//        s.d[i] = d[i] + t.d[i];
    s.d[ 0] = d[ 0] + t.d[ 0];    s.d[27] = d[27] + t.d[27];    s.d[54] = d[54] + t.d[54];
    s.d[ 1] = d[ 1] + t.d[ 1];    s.d[28] = d[28] + t.d[28];    s.d[55] = d[55] + t.d[55];
    s.d[ 2] = d[ 2] + t.d[ 2];    s.d[29] = d[29] + t.d[29];    s.d[56] = d[56] + t.d[56];
    s.d[ 3] = d[ 3] + t.d[ 3];    s.d[30] = d[30] + t.d[30];    s.d[57] = d[57] + t.d[57];
    s.d[ 4] = d[ 4] + t.d[ 4];    s.d[31] = d[31] + t.d[31];    s.d[58] = d[58] + t.d[58];
    s.d[ 5] = d[ 5] + t.d[ 5];    s.d[32] = d[32] + t.d[32];    s.d[59] = d[59] + t.d[59];
    s.d[ 6] = d[ 6] + t.d[ 6];    s.d[33] = d[33] + t.d[33];    s.d[60] = d[60] + t.d[60];
    s.d[ 7] = d[ 7] + t.d[ 7];    s.d[34] = d[34] + t.d[34];    s.d[61] = d[61] + t.d[61];
    s.d[ 8] = d[ 8] + t.d[ 8];    s.d[35] = d[35] + t.d[35];    s.d[62] = d[62] + t.d[62];
    s.d[ 9] = d[ 9] + t.d[ 9];    s.d[36] = d[36] + t.d[36];    s.d[63] = d[63] + t.d[63];
    s.d[10] = d[10] + t.d[10];    s.d[37] = d[37] + t.d[37];    s.d[64] = d[64] + t.d[64];
    s.d[11] = d[11] + t.d[11];    s.d[38] = d[38] + t.d[38];    s.d[65] = d[65] + t.d[65];
    s.d[12] = d[12] + t.d[12];    s.d[39] = d[39] + t.d[39];    s.d[66] = d[66] + t.d[66];
    s.d[13] = d[13] + t.d[13];    s.d[40] = d[40] + t.d[40];    s.d[67] = d[67] + t.d[67];
    s.d[14] = d[14] + t.d[14];    s.d[41] = d[41] + t.d[41];    s.d[68] = d[68] + t.d[68];
    s.d[15] = d[15] + t.d[15];    s.d[42] = d[42] + t.d[42];    s.d[69] = d[69] + t.d[69];
    s.d[16] = d[16] + t.d[16];    s.d[43] = d[43] + t.d[43];    s.d[70] = d[70] + t.d[70];
    s.d[17] = d[17] + t.d[17];    s.d[44] = d[44] + t.d[44];    s.d[71] = d[71] + t.d[71];
    s.d[18] = d[18] + t.d[18];    s.d[45] = d[45] + t.d[45];    s.d[72] = d[72] + t.d[72];
    s.d[19] = d[19] + t.d[19];    s.d[46] = d[46] + t.d[46];    s.d[73] = d[73] + t.d[73];
    s.d[20] = d[20] + t.d[20];    s.d[47] = d[47] + t.d[47];    s.d[74] = d[74] + t.d[74];
    s.d[21] = d[21] + t.d[21];    s.d[48] = d[48] + t.d[48];    s.d[75] = d[75] + t.d[75];
    s.d[22] = d[22] + t.d[22];    s.d[49] = d[49] + t.d[49];    s.d[76] = d[76] + t.d[76];
    s.d[23] = d[23] + t.d[23];    s.d[50] = d[50] + t.d[50];    s.d[77] = d[77] + t.d[77];
    s.d[24] = d[24] + t.d[24];    s.d[51] = d[51] + t.d[51];    s.d[78] = d[78] + t.d[78];
    s.d[25] = d[25] + t.d[25];    s.d[52] = d[52] + t.d[52];    s.d[79] = d[79] + t.d[79];
    s.d[26] = d[26] + t.d[26];    s.d[53] = d[53] + t.d[53];    s.d[80] = d[80] + t.d[80];
    return s;
}

//-----------------------------------------------------------------------------
// operator -
inline tens4d tens4d::operator - (const tens4d& t) const
{
    tens4d s;
//    for (int i=0; i<NNZ; i++)
//        s.d[i] = d[i] - t.d[i];
    s.d[ 0] = d[ 0] - t.d[ 0];    s.d[27] = d[27] - t.d[27];    s.d[54] = d[54] - t.d[54];
    s.d[ 1] = d[ 1] - t.d[ 1];    s.d[28] = d[28] - t.d[28];    s.d[55] = d[55] - t.d[55];
    s.d[ 2] = d[ 2] - t.d[ 2];    s.d[29] = d[29] - t.d[29];    s.d[56] = d[56] - t.d[56];
    s.d[ 3] = d[ 3] - t.d[ 3];    s.d[30] = d[30] - t.d[30];    s.d[57] = d[57] - t.d[57];
    s.d[ 4] = d[ 4] - t.d[ 4];    s.d[31] = d[31] - t.d[31];    s.d[58] = d[58] - t.d[58];
    s.d[ 5] = d[ 5] - t.d[ 5];    s.d[32] = d[32] - t.d[32];    s.d[59] = d[59] - t.d[59];
    s.d[ 6] = d[ 6] - t.d[ 6];    s.d[33] = d[33] - t.d[33];    s.d[60] = d[60] - t.d[60];
    s.d[ 7] = d[ 7] - t.d[ 7];    s.d[34] = d[34] - t.d[34];    s.d[61] = d[61] - t.d[61];
    s.d[ 8] = d[ 8] - t.d[ 8];    s.d[35] = d[35] - t.d[35];    s.d[62] = d[62] - t.d[62];
    s.d[ 9] = d[ 9] - t.d[ 9];    s.d[36] = d[36] - t.d[36];    s.d[63] = d[63] - t.d[63];
    s.d[10] = d[10] - t.d[10];    s.d[37] = d[37] - t.d[37];    s.d[64] = d[64] - t.d[64];
    s.d[11] = d[11] - t.d[11];    s.d[38] = d[38] - t.d[38];    s.d[65] = d[65] - t.d[65];
    s.d[12] = d[12] - t.d[12];    s.d[39] = d[39] - t.d[39];    s.d[66] = d[66] - t.d[66];
    s.d[13] = d[13] - t.d[13];    s.d[40] = d[40] - t.d[40];    s.d[67] = d[67] - t.d[67];
    s.d[14] = d[14] - t.d[14];    s.d[41] = d[41] - t.d[41];    s.d[68] = d[68] - t.d[68];
    s.d[15] = d[15] - t.d[15];    s.d[42] = d[42] - t.d[42];    s.d[69] = d[69] - t.d[69];
    s.d[16] = d[16] - t.d[16];    s.d[43] = d[43] - t.d[43];    s.d[70] = d[70] - t.d[70];
    s.d[17] = d[17] - t.d[17];    s.d[44] = d[44] - t.d[44];    s.d[71] = d[71] - t.d[71];
    s.d[18] = d[18] - t.d[18];    s.d[45] = d[45] - t.d[45];    s.d[72] = d[72] - t.d[72];
    s.d[19] = d[19] - t.d[19];    s.d[46] = d[46] - t.d[46];    s.d[73] = d[73] - t.d[73];
    s.d[20] = d[20] - t.d[20];    s.d[47] = d[47] - t.d[47];    s.d[74] = d[74] - t.d[74];
    s.d[21] = d[21] - t.d[21];    s.d[48] = d[48] - t.d[48];    s.d[75] = d[75] - t.d[75];
    s.d[22] = d[22] - t.d[22];    s.d[49] = d[49] - t.d[49];    s.d[76] = d[76] - t.d[76];
    s.d[23] = d[23] - t.d[23];    s.d[50] = d[50] - t.d[50];    s.d[77] = d[77] - t.d[77];
    s.d[24] = d[24] - t.d[24];    s.d[51] = d[51] - t.d[51];    s.d[78] = d[78] - t.d[78];
    s.d[25] = d[25] - t.d[25];    s.d[52] = d[52] - t.d[52];    s.d[79] = d[79] - t.d[79];
    s.d[26] = d[26] - t.d[26];    s.d[53] = d[53] - t.d[53];    s.d[80] = d[80] - t.d[80];
    return s;
}

//-----------------------------------------------------------------------------
// operator + tens4ds
inline tens4d tens4d::operator + (const tens4ds& t) const
{
    tens4d s;
    s.d[ 0] = d[ 0] + t.d[ 0];    s.d[27] = d[27] + t.d[ 6];    s.d[54] = d[54] + t.d[ 6];
    s.d[ 1] = d[ 1] + t.d[ 1];    s.d[28] = d[28] + t.d[ 7];    s.d[55] = d[55] + t.d[ 7];
    s.d[ 2] = d[ 2] + t.d[ 3];    s.d[29] = d[29] + t.d[ 8];    s.d[56] = d[56] + t.d[ 8];
    s.d[ 3] = d[ 3] + t.d[ 6];    s.d[30] = d[30] + t.d[ 9];    s.d[57] = d[57] + t.d[ 9];
    s.d[ 4] = d[ 4] + t.d[10];    s.d[31] = d[31] + t.d[13];    s.d[58] = d[58] + t.d[13];
    s.d[ 5] = d[ 5] + t.d[15];    s.d[32] = d[32] + t.d[18];    s.d[59] = d[59] + t.d[18];
    s.d[ 6] = d[ 6] + t.d[ 6];    s.d[33] = d[33] + t.d[ 9];    s.d[60] = d[60] + t.d[9];
    s.d[ 7] = d[ 7] + t.d[10];    s.d[34] = d[34] + t.d[13];    s.d[61] = d[61] + t.d[13];
    s.d[ 8] = d[ 8] + t.d[15];    s.d[35] = d[35] + t.d[18];    s.d[62] = d[62] + t.d[18];
    s.d[ 9] = d[ 9] + t.d[ 1];    s.d[36] = d[36] + t.d[10];    s.d[63] = d[63] + t.d[10];
    s.d[10] = d[10] + t.d[ 2];    s.d[37] = d[37] + t.d[11];    s.d[64] = d[64] + t.d[11];
    s.d[11] = d[11] + t.d[ 4];    s.d[38] = d[38] + t.d[12];    s.d[65] = d[65] + t.d[12];
    s.d[12] = d[12] + t.d[ 7];    s.d[39] = d[39] + t.d[13];    s.d[66] = d[66] + t.d[13];
    s.d[13] = d[13] + t.d[11];    s.d[40] = d[40] + t.d[14];    s.d[67] = d[67] + t.d[14];
    s.d[14] = d[14] + t.d[16];    s.d[41] = d[41] + t.d[19];    s.d[68] = d[68] + t.d[19];
    s.d[15] = d[15] + t.d[ 7];    s.d[42] = d[42] + t.d[13];    s.d[69] = d[69] + t.d[13];
    s.d[16] = d[16] + t.d[11];    s.d[43] = d[43] + t.d[14];    s.d[70] = d[70] + t.d[14];
    s.d[17] = d[17] + t.d[16];    s.d[44] = d[44] + t.d[19];    s.d[71] = d[71] + t.d[19];
    s.d[18] = d[18] + t.d[ 3];    s.d[45] = d[45] + t.d[15];    s.d[72] = d[72] + t.d[15];
    s.d[19] = d[19] + t.d[ 4];    s.d[46] = d[46] + t.d[16];    s.d[73] = d[73] + t.d[16];
    s.d[20] = d[20] + t.d[ 5];    s.d[47] = d[47] + t.d[17];    s.d[74] = d[74] + t.d[17];
    s.d[21] = d[21] + t.d[ 8];    s.d[48] = d[48] + t.d[18];    s.d[75] = d[75] + t.d[18];
    s.d[22] = d[22] + t.d[12];    s.d[49] = d[49] + t.d[19];    s.d[76] = d[76] + t.d[19];
    s.d[23] = d[23] + t.d[17];    s.d[50] = d[50] + t.d[20];    s.d[77] = d[77] + t.d[20];
    s.d[24] = d[24] + t.d[ 8];    s.d[51] = d[51] + t.d[18];    s.d[78] = d[78] + t.d[18];
    s.d[25] = d[25] + t.d[12];    s.d[52] = d[52] + t.d[19];    s.d[79] = d[79] + t.d[19];
    s.d[26] = d[26] + t.d[17];    s.d[53] = d[53] + t.d[20];    s.d[80] = d[80] + t.d[20];
    return s;
}

//-----------------------------------------------------------------------------
// operator - tens4ds
inline tens4d tens4d::operator - (const tens4ds& t) const
{
    tens4d s;
    s.d[ 0] = d[ 0] - t.d[ 0];    s.d[27] = d[27] - t.d[ 6];    s.d[54] = d[54] - t.d[ 6];
    s.d[ 1] = d[ 1] - t.d[ 1];    s.d[28] = d[28] - t.d[ 7];    s.d[55] = d[55] - t.d[ 7];
    s.d[ 2] = d[ 2] - t.d[ 3];    s.d[29] = d[29] - t.d[ 8];    s.d[56] = d[56] - t.d[ 8];
    s.d[ 3] = d[ 3] - t.d[ 6];    s.d[30] = d[30] - t.d[ 9];    s.d[57] = d[57] - t.d[ 9];
    s.d[ 4] = d[ 4] - t.d[10];    s.d[31] = d[31] - t.d[13];    s.d[58] = d[58] - t.d[13];
    s.d[ 5] = d[ 5] - t.d[15];    s.d[32] = d[32] - t.d[18];    s.d[59] = d[59] - t.d[18];
    s.d[ 6] = d[ 6] - t.d[ 6];    s.d[33] = d[33] - t.d[ 9];    s.d[60] = d[60] - t.d[9];
    s.d[ 7] = d[ 7] - t.d[10];    s.d[34] = d[34] - t.d[13];    s.d[61] = d[61] - t.d[13];
    s.d[ 8] = d[ 8] - t.d[15];    s.d[35] = d[35] - t.d[18];    s.d[62] = d[62] - t.d[18];
    s.d[ 9] = d[ 9] - t.d[ 1];    s.d[36] = d[36] - t.d[10];    s.d[63] = d[63] - t.d[10];
    s.d[10] = d[10] - t.d[ 2];    s.d[37] = d[37] - t.d[11];    s.d[64] = d[64] - t.d[11];
    s.d[11] = d[11] - t.d[ 4];    s.d[38] = d[38] - t.d[12];    s.d[65] = d[65] - t.d[12];
    s.d[12] = d[12] - t.d[ 7];    s.d[39] = d[39] - t.d[13];    s.d[66] = d[66] - t.d[13];
    s.d[13] = d[13] - t.d[11];    s.d[40] = d[40] - t.d[14];    s.d[67] = d[67] - t.d[14];
    s.d[14] = d[14] - t.d[16];    s.d[41] = d[41] - t.d[19];    s.d[68] = d[68] - t.d[19];
    s.d[15] = d[15] - t.d[ 7];    s.d[42] = d[42] - t.d[13];    s.d[69] = d[69] - t.d[13];
    s.d[16] = d[16] - t.d[11];    s.d[43] = d[43] - t.d[14];    s.d[70] = d[70] - t.d[14];
    s.d[17] = d[17] - t.d[16];    s.d[44] = d[44] - t.d[19];    s.d[71] = d[71] - t.d[19];
    s.d[18] = d[18] - t.d[ 3];    s.d[45] = d[45] - t.d[15];    s.d[72] = d[72] - t.d[15];
    s.d[19] = d[19] - t.d[ 4];    s.d[46] = d[46] - t.d[16];    s.d[73] = d[73] - t.d[16];
    s.d[20] = d[20] - t.d[ 5];    s.d[47] = d[47] - t.d[17];    s.d[74] = d[74] - t.d[17];
    s.d[21] = d[21] - t.d[ 8];    s.d[48] = d[48] - t.d[18];    s.d[75] = d[75] - t.d[18];
    s.d[22] = d[22] - t.d[12];    s.d[49] = d[49] - t.d[19];    s.d[76] = d[76] - t.d[19];
    s.d[23] = d[23] - t.d[17];    s.d[50] = d[50] - t.d[20];    s.d[77] = d[77] - t.d[20];
    s.d[24] = d[24] - t.d[ 8];    s.d[51] = d[51] - t.d[18];    s.d[78] = d[78] - t.d[18];
    s.d[25] = d[25] - t.d[12];    s.d[52] = d[52] - t.d[19];    s.d[79] = d[79] - t.d[19];
    s.d[26] = d[26] - t.d[17];    s.d[53] = d[53] - t.d[20];    s.d[80] = d[80] - t.d[20];
    return s;
}

//-----------------------------------------------------------------------------
// operator *
inline tens4d tens4d::operator * (double g) const
{
    tens4d s;
//    for (int i=0; i<NNZ; i++)
//        s.d[i] = g*d[i];
    s.d[ 0] = g*d[ 0];    s.d[27] = g*d[27];    s.d[54] = g*d[54];
    s.d[ 1] = g*d[ 1];    s.d[28] = g*d[28];    s.d[55] = g*d[55];
    s.d[ 2] = g*d[ 2];    s.d[29] = g*d[29];    s.d[56] = g*d[56];
    s.d[ 3] = g*d[ 3];    s.d[30] = g*d[30];    s.d[57] = g*d[57];
    s.d[ 4] = g*d[ 4];    s.d[31] = g*d[31];    s.d[58] = g*d[58];
    s.d[ 5] = g*d[ 5];    s.d[32] = g*d[32];    s.d[59] = g*d[59];
    s.d[ 6] = g*d[ 6];    s.d[33] = g*d[33];    s.d[60] = g*d[60];
    s.d[ 7] = g*d[ 7];    s.d[34] = g*d[34];    s.d[61] = g*d[61];
    s.d[ 8] = g*d[ 8];    s.d[35] = g*d[35];    s.d[62] = g*d[62];
    s.d[ 9] = g*d[ 9];    s.d[36] = g*d[36];    s.d[63] = g*d[63];
    s.d[10] = g*d[10];    s.d[37] = g*d[37];    s.d[64] = g*d[64];
    s.d[11] = g*d[11];    s.d[38] = g*d[38];    s.d[65] = g*d[65];
    s.d[12] = g*d[12];    s.d[39] = g*d[39];    s.d[66] = g*d[66];
    s.d[13] = g*d[13];    s.d[40] = g*d[40];    s.d[67] = g*d[67];
    s.d[14] = g*d[14];    s.d[41] = g*d[41];    s.d[68] = g*d[68];
    s.d[15] = g*d[15];    s.d[42] = g*d[42];    s.d[69] = g*d[69];
    s.d[16] = g*d[16];    s.d[43] = g*d[43];    s.d[70] = g*d[70];
    s.d[17] = g*d[17];    s.d[44] = g*d[44];    s.d[71] = g*d[71];
    s.d[18] = g*d[18];    s.d[45] = g*d[45];    s.d[72] = g*d[72];
    s.d[19] = g*d[19];    s.d[46] = g*d[46];    s.d[73] = g*d[73];
    s.d[20] = g*d[20];    s.d[47] = g*d[47];    s.d[74] = g*d[74];
    s.d[21] = g*d[21];    s.d[48] = g*d[48];    s.d[75] = g*d[75];
    s.d[22] = g*d[22];    s.d[49] = g*d[49];    s.d[76] = g*d[76];
    s.d[23] = g*d[23];    s.d[50] = g*d[50];    s.d[77] = g*d[77];
    s.d[24] = g*d[24];    s.d[51] = g*d[51];    s.d[78] = g*d[78];
    s.d[25] = g*d[25];    s.d[52] = g*d[52];    s.d[79] = g*d[79];
    s.d[26] = g*d[26];    s.d[53] = g*d[53];    s.d[80] = g*d[80];
    return s;
}

//-----------------------------------------------------------------------------
// operator /
inline tens4d tens4d::operator / (double g) const
{
    tens4d s;
//    for (int i=0; i<NNZ; i++)
//        s.d[i] = d[i]/g;
    s.d[ 0] = d[ 0]/g;    s.d[27] = d[27]/g;    s.d[54] = d[54]/g;
    s.d[ 1] = d[ 1]/g;    s.d[28] = d[28]/g;    s.d[55] = d[55]/g;
    s.d[ 2] = d[ 2]/g;    s.d[29] = d[29]/g;    s.d[56] = d[56]/g;
    s.d[ 3] = d[ 3]/g;    s.d[30] = d[30]/g;    s.d[57] = d[57]/g;
    s.d[ 4] = d[ 4]/g;    s.d[31] = d[31]/g;    s.d[58] = d[58]/g;
    s.d[ 5] = d[ 5]/g;    s.d[32] = d[32]/g;    s.d[59] = d[59]/g;
    s.d[ 6] = d[ 6]/g;    s.d[33] = d[33]/g;    s.d[60] = d[60]/g;
    s.d[ 7] = d[ 7]/g;    s.d[34] = d[34]/g;    s.d[61] = d[61]/g;
    s.d[ 8] = d[ 8]/g;    s.d[35] = d[35]/g;    s.d[62] = d[62]/g;
    s.d[ 9] = d[ 9]/g;    s.d[36] = d[36]/g;    s.d[63] = d[63]/g;
    s.d[10] = d[10]/g;    s.d[37] = d[37]/g;    s.d[64] = d[64]/g;
    s.d[11] = d[11]/g;    s.d[38] = d[38]/g;    s.d[65] = d[65]/g;
    s.d[12] = d[12]/g;    s.d[39] = d[39]/g;    s.d[66] = d[66]/g;
    s.d[13] = d[13]/g;    s.d[40] = d[40]/g;    s.d[67] = d[67]/g;
    s.d[14] = d[14]/g;    s.d[41] = d[41]/g;    s.d[68] = d[68]/g;
    s.d[15] = d[15]/g;    s.d[42] = d[42]/g;    s.d[69] = d[69]/g;
    s.d[16] = d[16]/g;    s.d[43] = d[43]/g;    s.d[70] = d[70]/g;
    s.d[17] = d[17]/g;    s.d[44] = d[44]/g;    s.d[71] = d[71]/g;
    s.d[18] = d[18]/g;    s.d[45] = d[45]/g;    s.d[72] = d[72]/g;
    s.d[19] = d[19]/g;    s.d[46] = d[46]/g;    s.d[73] = d[73]/g;
    s.d[20] = d[20]/g;    s.d[47] = d[47]/g;    s.d[74] = d[74]/g;
    s.d[21] = d[21]/g;    s.d[48] = d[48]/g;    s.d[75] = d[75]/g;
    s.d[22] = d[22]/g;    s.d[49] = d[49]/g;    s.d[76] = d[76]/g;
    s.d[23] = d[23]/g;    s.d[50] = d[50]/g;    s.d[77] = d[77]/g;
    s.d[24] = d[24]/g;    s.d[51] = d[51]/g;    s.d[78] = d[78]/g;
    s.d[25] = d[25]/g;    s.d[52] = d[52]/g;    s.d[79] = d[79]/g;
    s.d[26] = d[26]/g;    s.d[53] = d[53]/g;    s.d[80] = d[80]/g;
    return s;
}

//-----------------------------------------------------------------------------
// assignment operator +=
inline tens4d& tens4d::operator += (const tens4d& t)
{
//    for (int i=0; i<NNZ; i++)
//        d[i] += t.d[i];
    d[ 0] += t.d[ 0];    d[27] += t.d[27];    d[54] += t.d[54];
    d[ 1] += t.d[ 1];    d[28] += t.d[28];    d[55] += t.d[55];
    d[ 2] += t.d[ 2];    d[29] += t.d[29];    d[56] += t.d[56];
    d[ 3] += t.d[ 3];    d[30] += t.d[30];    d[57] += t.d[57];
    d[ 4] += t.d[ 4];    d[31] += t.d[31];    d[58] += t.d[58];
    d[ 5] += t.d[ 5];    d[32] += t.d[32];    d[59] += t.d[59];
    d[ 6] += t.d[ 6];    d[33] += t.d[33];    d[60] += t.d[60];
    d[ 7] += t.d[ 7];    d[34] += t.d[34];    d[61] += t.d[61];
    d[ 8] += t.d[ 8];    d[35] += t.d[35];    d[62] += t.d[62];
    d[ 9] += t.d[ 9];    d[36] += t.d[36];    d[63] += t.d[63];
    d[10] += t.d[10];    d[37] += t.d[37];    d[64] += t.d[64];
    d[11] += t.d[11];    d[38] += t.d[38];    d[65] += t.d[65];
    d[12] += t.d[12];    d[39] += t.d[39];    d[66] += t.d[66];
    d[13] += t.d[13];    d[40] += t.d[40];    d[67] += t.d[67];
    d[14] += t.d[14];    d[41] += t.d[41];    d[68] += t.d[68];
    d[15] += t.d[15];    d[42] += t.d[42];    d[69] += t.d[69];
    d[16] += t.d[16];    d[43] += t.d[43];    d[70] += t.d[70];
    d[17] += t.d[17];    d[44] += t.d[44];    d[71] += t.d[71];
    d[18] += t.d[18];    d[45] += t.d[45];    d[72] += t.d[72];
    d[19] += t.d[19];    d[46] += t.d[46];    d[73] += t.d[73];
    d[20] += t.d[20];    d[47] += t.d[47];    d[74] += t.d[74];
    d[21] += t.d[21];    d[48] += t.d[48];    d[75] += t.d[75];
    d[22] += t.d[22];    d[49] += t.d[49];    d[76] += t.d[76];
    d[23] += t.d[23];    d[50] += t.d[50];    d[77] += t.d[77];
    d[24] += t.d[24];    d[51] += t.d[51];    d[78] += t.d[78];
    d[25] += t.d[25];    d[52] += t.d[52];    d[79] += t.d[79];
    d[26] += t.d[26];    d[53] += t.d[53];    d[80] += t.d[80];
    return (*this);
}

//-----------------------------------------------------------------------------
// assignment operator -=
inline tens4d& tens4d::operator -= (const tens4d& t)
{
//    for (int i=0; i<NNZ; i++)
//        d[i] -= t.d[i];
    d[ 0] -= t.d[ 0];    d[27] -= t.d[27];    d[54] -= t.d[54];
    d[ 1] -= t.d[ 1];    d[28] -= t.d[28];    d[55] -= t.d[55];
    d[ 2] -= t.d[ 2];    d[29] -= t.d[29];    d[56] -= t.d[56];
    d[ 3] -= t.d[ 3];    d[30] -= t.d[30];    d[57] -= t.d[57];
    d[ 4] -= t.d[ 4];    d[31] -= t.d[31];    d[58] -= t.d[58];
    d[ 5] -= t.d[ 5];    d[32] -= t.d[32];    d[59] -= t.d[59];
    d[ 6] -= t.d[ 6];    d[33] -= t.d[33];    d[60] -= t.d[60];
    d[ 7] -= t.d[ 7];    d[34] -= t.d[34];    d[61] -= t.d[61];
    d[ 8] -= t.d[ 8];    d[35] -= t.d[35];    d[62] -= t.d[62];
    d[ 9] -= t.d[ 9];    d[36] -= t.d[36];    d[63] -= t.d[63];
    d[10] -= t.d[10];    d[37] -= t.d[37];    d[64] -= t.d[64];
    d[11] -= t.d[11];    d[38] -= t.d[38];    d[65] -= t.d[65];
    d[12] -= t.d[12];    d[39] -= t.d[39];    d[66] -= t.d[66];
    d[13] -= t.d[13];    d[40] -= t.d[40];    d[67] -= t.d[67];
    d[14] -= t.d[14];    d[41] -= t.d[41];    d[68] -= t.d[68];
    d[15] -= t.d[15];    d[42] -= t.d[42];    d[69] -= t.d[69];
    d[16] -= t.d[16];    d[43] -= t.d[43];    d[70] -= t.d[70];
    d[17] -= t.d[17];    d[44] -= t.d[44];    d[71] -= t.d[71];
    d[18] -= t.d[18];    d[45] -= t.d[45];    d[72] -= t.d[72];
    d[19] -= t.d[19];    d[46] -= t.d[46];    d[73] -= t.d[73];
    d[20] -= t.d[20];    d[47] -= t.d[47];    d[74] -= t.d[74];
    d[21] -= t.d[21];    d[48] -= t.d[48];    d[75] -= t.d[75];
    d[22] -= t.d[22];    d[49] -= t.d[49];    d[76] -= t.d[76];
    d[23] -= t.d[23];    d[50] -= t.d[50];    d[77] -= t.d[77];
    d[24] -= t.d[24];    d[51] -= t.d[51];    d[78] -= t.d[78];
    d[25] -= t.d[25];    d[52] -= t.d[52];    d[79] -= t.d[79];
    d[26] -= t.d[26];    d[53] -= t.d[53];    d[80] -= t.d[80];
    return (*this);
}

//-----------------------------------------------------------------------------
// assignment operator += tens4ds
inline tens4d& tens4d::operator += (const tens4ds& t)
{
    d[ 0] += t.d[ 0];    d[27] += t.d[ 6];    d[54] += t.d[ 6];
    d[ 1] += t.d[ 1];    d[28] += t.d[ 7];    d[55] += t.d[ 7];
    d[ 2] += t.d[ 3];    d[29] += t.d[ 8];    d[56] += t.d[ 8];
    d[ 3] += t.d[ 6];    d[30] += t.d[ 9];    d[57] += t.d[ 9];
    d[ 4] += t.d[10];    d[31] += t.d[13];    d[58] += t.d[13];
    d[ 5] += t.d[15];    d[32] += t.d[18];    d[59] += t.d[18];
    d[ 6] += t.d[ 6];    d[33] += t.d[ 9];    d[60] += t.d[ 9];
    d[ 7] += t.d[10];    d[34] += t.d[13];    d[61] += t.d[13];
    d[ 8] += t.d[15];    d[35] += t.d[18];    d[62] += t.d[18];
    d[ 9] += t.d[ 1];    d[36] += t.d[10];    d[63] += t.d[10];
    d[10] += t.d[ 2];    d[37] += t.d[11];    d[64] += t.d[11];
    d[11] += t.d[ 4];    d[38] += t.d[12];    d[65] += t.d[12];
    d[12] += t.d[ 7];    d[39] += t.d[13];    d[66] += t.d[13];
    d[13] += t.d[11];    d[40] += t.d[14];    d[67] += t.d[14];
    d[14] += t.d[16];    d[41] += t.d[19];    d[68] += t.d[19];
    d[15] += t.d[ 7];    d[42] += t.d[13];    d[69] += t.d[13];
    d[16] += t.d[11];    d[43] += t.d[14];    d[70] += t.d[14];
    d[17] += t.d[16];    d[44] += t.d[19];    d[71] += t.d[19];
    d[18] += t.d[ 3];    d[45] += t.d[15];    d[72] += t.d[15];
    d[19] += t.d[ 4];    d[46] += t.d[16];    d[73] += t.d[16];
    d[20] += t.d[ 5];    d[47] += t.d[17];    d[74] += t.d[17];
    d[21] += t.d[ 8];    d[48] += t.d[18];    d[75] += t.d[18];
    d[22] += t.d[12];    d[49] += t.d[19];    d[76] += t.d[19];
    d[23] += t.d[17];    d[50] += t.d[20];    d[77] += t.d[20];
    d[24] += t.d[ 8];    d[51] += t.d[18];    d[78] += t.d[18];
    d[25] += t.d[12];    d[52] += t.d[19];    d[79] += t.d[19];
    d[26] += t.d[17];    d[53] += t.d[20];    d[80] += t.d[20];
    return (*this);
}

//-----------------------------------------------------------------------------
// assignment operator -= tens4ds
inline tens4d& tens4d::operator -= (const tens4ds& t)
{
    d[ 0] -= t.d[ 0];    d[27] -= t.d[ 6];    d[54] -= t.d[ 6];
    d[ 1] -= t.d[ 1];    d[28] -= t.d[ 7];    d[55] -= t.d[ 7];
    d[ 2] -= t.d[ 3];    d[29] -= t.d[ 8];    d[56] -= t.d[ 8];
    d[ 3] -= t.d[ 6];    d[30] -= t.d[ 9];    d[57] -= t.d[ 9];
    d[ 4] -= t.d[10];    d[31] -= t.d[13];    d[58] -= t.d[13];
    d[ 5] -= t.d[15];    d[32] -= t.d[18];    d[59] -= t.d[18];
    d[ 6] -= t.d[ 6];    d[33] -= t.d[ 9];    d[60] -= t.d[ 9];
    d[ 7] -= t.d[10];    d[34] -= t.d[13];    d[61] -= t.d[13];
    d[ 8] -= t.d[15];    d[35] -= t.d[18];    d[62] -= t.d[18];
    d[ 9] -= t.d[ 1];    d[36] -= t.d[10];    d[63] -= t.d[10];
    d[10] -= t.d[ 2];    d[37] -= t.d[11];    d[64] -= t.d[11];
    d[11] -= t.d[ 4];    d[38] -= t.d[12];    d[65] -= t.d[12];
    d[12] -= t.d[ 7];    d[39] -= t.d[13];    d[66] -= t.d[13];
    d[13] -= t.d[11];    d[40] -= t.d[14];    d[67] -= t.d[14];
    d[14] -= t.d[16];    d[41] -= t.d[19];    d[68] -= t.d[19];
    d[15] -= t.d[ 7];    d[42] -= t.d[13];    d[69] -= t.d[13];
    d[16] -= t.d[11];    d[43] -= t.d[14];    d[70] -= t.d[14];
    d[17] -= t.d[16];    d[44] -= t.d[19];    d[71] -= t.d[19];
    d[18] -= t.d[ 3];    d[45] -= t.d[15];    d[72] -= t.d[15];
    d[19] -= t.d[ 4];    d[46] -= t.d[16];    d[73] -= t.d[16];
    d[20] -= t.d[ 5];    d[47] -= t.d[17];    d[74] -= t.d[17];
    d[21] -= t.d[ 8];    d[48] -= t.d[18];    d[75] -= t.d[18];
    d[22] -= t.d[12];    d[49] -= t.d[19];    d[76] -= t.d[19];
    d[23] -= t.d[17];    d[50] -= t.d[20];    d[77] -= t.d[20];
    d[24] -= t.d[ 8];    d[51] -= t.d[18];    d[78] -= t.d[18];
    d[25] -= t.d[12];    d[52] -= t.d[19];    d[79] -= t.d[19];
    d[26] -= t.d[17];    d[53] -= t.d[20];    d[80] -= t.d[20];
    return (*this);
}

//-----------------------------------------------------------------------------
// assignment operator *=
inline tens4d& tens4d::operator *= (double g)
{
//    for (int i=0; i<NNZ; i++)
//        d[i] *= g;
    d[ 0] *= g;    d[27] *= g;    d[54] *= g;
    d[ 1] *= g;    d[28] *= g;    d[55] *= g;
    d[ 2] *= g;    d[29] *= g;    d[56] *= g;
    d[ 3] *= g;    d[30] *= g;    d[57] *= g;
    d[ 4] *= g;    d[31] *= g;    d[58] *= g;
    d[ 5] *= g;    d[32] *= g;    d[59] *= g;
    d[ 6] *= g;    d[33] *= g;    d[60] *= g;
    d[ 7] *= g;    d[34] *= g;    d[61] *= g;
    d[ 8] *= g;    d[35] *= g;    d[62] *= g;
    d[ 9] *= g;    d[36] *= g;    d[63] *= g;
    d[10] *= g;    d[37] *= g;    d[64] *= g;
    d[11] *= g;    d[38] *= g;    d[65] *= g;
    d[12] *= g;    d[39] *= g;    d[66] *= g;
    d[13] *= g;    d[40] *= g;    d[67] *= g;
    d[14] *= g;    d[41] *= g;    d[68] *= g;
    d[15] *= g;    d[42] *= g;    d[69] *= g;
    d[16] *= g;    d[43] *= g;    d[70] *= g;
    d[17] *= g;    d[44] *= g;    d[71] *= g;
    d[18] *= g;    d[45] *= g;    d[72] *= g;
    d[19] *= g;    d[46] *= g;    d[73] *= g;
    d[20] *= g;    d[47] *= g;    d[74] *= g;
    d[21] *= g;    d[48] *= g;    d[75] *= g;
    d[22] *= g;    d[49] *= g;    d[76] *= g;
    d[23] *= g;    d[50] *= g;    d[77] *= g;
    d[24] *= g;    d[51] *= g;    d[78] *= g;
    d[25] *= g;    d[52] *= g;    d[79] *= g;
    d[26] *= g;    d[53] *= g;    d[80] *= g;
    return (*this);
}

//-----------------------------------------------------------------------------
// assignment operator /=
inline tens4d& tens4d::operator /= (double g)
{
//    for (int i=0; i<NNZ; i++)
//        d[i] /= g;
    d[ 0] /= g;    d[27] /= g;    d[54] /= g;
    d[ 1] /= g;    d[28] /= g;    d[55] /= g;
    d[ 2] /= g;    d[29] /= g;    d[56] /= g;
    d[ 3] /= g;    d[30] /= g;    d[57] /= g;
    d[ 4] /= g;    d[31] /= g;    d[58] /= g;
    d[ 5] /= g;    d[32] /= g;    d[59] /= g;
    d[ 6] /= g;    d[33] /= g;    d[60] /= g;
    d[ 7] /= g;    d[34] /= g;    d[61] /= g;
    d[ 8] /= g;    d[35] /= g;    d[62] /= g;
    d[ 9] /= g;    d[36] /= g;    d[63] /= g;
    d[10] /= g;    d[37] /= g;    d[64] /= g;
    d[11] /= g;    d[38] /= g;    d[65] /= g;
    d[12] /= g;    d[39] /= g;    d[66] /= g;
    d[13] /= g;    d[40] /= g;    d[67] /= g;
    d[14] /= g;    d[41] /= g;    d[68] /= g;
    d[15] /= g;    d[42] /= g;    d[69] /= g;
    d[16] /= g;    d[43] /= g;    d[70] /= g;
    d[17] /= g;    d[44] /= g;    d[71] /= g;
    d[18] /= g;    d[45] /= g;    d[72] /= g;
    d[19] /= g;    d[46] /= g;    d[73] /= g;
    d[20] /= g;    d[47] /= g;    d[74] /= g;
    d[21] /= g;    d[48] /= g;    d[75] /= g;
    d[22] /= g;    d[49] /= g;    d[76] /= g;
    d[23] /= g;    d[50] /= g;    d[77] /= g;
    d[24] /= g;    d[51] /= g;    d[78] /= g;
    d[25] /= g;    d[52] /= g;    d[79] /= g;
    d[26] /= g;    d[53] /= g;    d[80] /= g;
    return (*this);
}

//-----------------------------------------------------------------------------
// unary operator -
inline tens4d tens4d::operator - () const
{
    tens4d s;
    s.d[ 0] = -d[ 0];    s.d[27] = -d[27];    s.d[54] = -d[54];
    s.d[ 1] = -d[ 1];    s.d[28] = -d[28];    s.d[55] = -d[55];
    s.d[ 2] = -d[ 2];    s.d[29] = -d[29];    s.d[56] = -d[56];
    s.d[ 3] = -d[ 3];    s.d[30] = -d[30];    s.d[57] = -d[57];
    s.d[ 4] = -d[ 4];    s.d[31] = -d[31];    s.d[58] = -d[58];
    s.d[ 5] = -d[ 5];    s.d[32] = -d[32];    s.d[59] = -d[59];
    s.d[ 6] = -d[ 6];    s.d[33] = -d[33];    s.d[60] = -d[60];
    s.d[ 7] = -d[ 7];    s.d[34] = -d[34];    s.d[61] = -d[61];
    s.d[ 8] = -d[ 8];    s.d[35] = -d[35];    s.d[62] = -d[62];
    s.d[ 9] = -d[ 9];    s.d[36] = -d[36];    s.d[63] = -d[63];
    s.d[10] = -d[10];    s.d[37] = -d[37];    s.d[64] = -d[64];
    s.d[11] = -d[11];    s.d[38] = -d[38];    s.d[65] = -d[65];
    s.d[12] = -d[12];    s.d[39] = -d[39];    s.d[66] = -d[66];
    s.d[13] = -d[13];    s.d[40] = -d[40];    s.d[67] = -d[67];
    s.d[14] = -d[14];    s.d[41] = -d[41];    s.d[68] = -d[68];
    s.d[15] = -d[15];    s.d[42] = -d[42];    s.d[69] = -d[69];
    s.d[16] = -d[16];    s.d[43] = -d[43];    s.d[70] = -d[70];
    s.d[17] = -d[17];    s.d[44] = -d[44];    s.d[71] = -d[71];
    s.d[18] = -d[18];    s.d[45] = -d[45];    s.d[72] = -d[72];
    s.d[19] = -d[19];    s.d[46] = -d[46];    s.d[73] = -d[73];
    s.d[20] = -d[20];    s.d[47] = -d[47];    s.d[74] = -d[74];
    s.d[21] = -d[21];    s.d[48] = -d[48];    s.d[75] = -d[75];
    s.d[22] = -d[22];    s.d[49] = -d[49];    s.d[76] = -d[76];
    s.d[23] = -d[23];    s.d[50] = -d[50];    s.d[77] = -d[77];
    s.d[24] = -d[24];    s.d[51] = -d[51];    s.d[78] = -d[78];
    s.d[25] = -d[25];    s.d[52] = -d[52];    s.d[79] = -d[79];
    s.d[26] = -d[26];    s.d[53] = -d[53];    s.d[80] = -d[80];
    return s;
}

//-----------------------------------------------------------------------------
// double contraction of general 4th-order tensor with a general 2nd-order tensor
// Aij = Dijkl Mkl
inline mat3d tens4d::dot(const mat3d &m) const
{
    mat3d a;
    
    a(0,0) = d[0]*m(0,0) + d[27]*m(0,1) + d[72]*m(0,2) + d[54]*m(1,0) + d[ 9]*m(1,1) + d[36]*m(1,2) + d[45]*m(2,0) + d[63]*m(2,1) + d[18]*m(2,2);
    a(0,1) = d[1]*m(0,0) + d[28]*m(0,1) + d[73]*m(0,2) + d[55]*m(1,0) + d[10]*m(1,1) + d[37]*m(1,2) + d[46]*m(2,0) + d[64]*m(2,1) + d[19]*m(2,2);
    a(0,2) = d[2]*m(0,0) + d[29]*m(0,1) + d[74]*m(0,2) + d[56]*m(1,0) + d[11]*m(1,1) + d[38]*m(1,2) + d[47]*m(2,0) + d[65]*m(2,1) + d[20]*m(2,2);
    a(1,0) = d[3]*m(0,0) + d[30]*m(0,1) + d[75]*m(0,2) + d[57]*m(1,0) + d[12]*m(1,1) + d[39]*m(1,2) + d[48]*m(2,0) + d[66]*m(2,1) + d[21]*m(2,2);
    a(1,1) = d[4]*m(0,0) + d[31]*m(0,1) + d[76]*m(0,2) + d[58]*m(1,0) + d[13]*m(1,1) + d[40]*m(1,2) + d[49]*m(2,0) + d[67]*m(2,1) + d[22]*m(2,2);
    a(1,2) = d[5]*m(0,0) + d[32]*m(0,1) + d[77]*m(0,2) + d[59]*m(1,0) + d[14]*m(1,1) + d[41]*m(1,2) + d[50]*m(2,0) + d[68]*m(2,1) + d[23]*m(2,2);
    a(2,0) = d[6]*m(0,0) + d[33]*m(0,1) + d[78]*m(0,2) + d[60]*m(1,0) + d[15]*m(1,1) + d[42]*m(1,2) + d[51]*m(2,0) + d[69]*m(2,1) + d[24]*m(2,2);
    a(2,1) = d[7]*m(0,0) + d[34]*m(0,1) + d[79]*m(0,2) + d[61]*m(1,0) + d[16]*m(1,1) + d[43]*m(1,2) + d[52]*m(2,0) + d[70]*m(2,1) + d[25]*m(2,2);
    a(2,2) = d[8]*m(0,0) + d[35]*m(0,1) + d[80]*m(0,2) + d[62]*m(1,0) + d[17]*m(1,1) + d[44]*m(1,2) + d[53]*m(2,0) + d[71]*m(2,1) + d[26]*m(2,2);
    
    return a;
}

// single dot product with 2nd order tensor
inline tens4d tens4d::sdot(const mat3d& m) const
{
    tens4d s;
    s.d[ 0] = d[0]*m(0,0) + d[27]*m(1,0) + d[72]*m(2,0);
    s.d[ 1] = d[1]*m(0,0) + d[28]*m(1,0) + d[73]*m(2,0);
    s.d[ 2] = d[2]*m(0,0) + d[29]*m(1,0) + d[74]*m(2,0);
    s.d[ 3] = d[3]*m(0,0) + d[30]*m(1,0) + d[75]*m(2,0);
    s.d[ 4] = d[4]*m(0,0) + d[31]*m(1,0) + d[76]*m(2,0);
    s.d[ 5] = d[5]*m(0,0) + d[32]*m(1,0) + d[77]*m(2,0);
    s.d[ 6] = d[6]*m(0,0) + d[33]*m(1,0) + d[78]*m(2,0);
    s.d[ 7] = d[7]*m(0,0) + d[34]*m(1,0) + d[79]*m(2,0);
    s.d[ 8] = d[8]*m(0,0) + d[35]*m(1,0) + d[80]*m(2,0);

    s.d[ 9] = d[54]*m(0,1) + d[ 9]*m(1,1) + d[36]*m(2,1);
    s.d[10] = d[55]*m(0,1) + d[10]*m(1,1) + d[37]*m(2,1);
    s.d[11] = d[56]*m(0,1) + d[11]*m(1,1) + d[38]*m(2,1);
    s.d[12] = d[57]*m(0,1) + d[12]*m(1,1) + d[39]*m(2,1);
    s.d[13] = d[58]*m(0,1) + d[13]*m(1,1) + d[40]*m(2,1);
    s.d[14] = d[59]*m(0,1) + d[14]*m(1,1) + d[41]*m(2,1);
    s.d[15] = d[60]*m(0,1) + d[15]*m(1,1) + d[42]*m(2,1);
    s.d[16] = d[61]*m(0,1) + d[16]*m(1,1) + d[43]*m(2,1);
    s.d[17] = d[62]*m(0,1) + d[17]*m(1,1) + d[44]*m(2,1);

    s.d[18] = d[45]*m(0,2) + d[63]*m(1,2) + d[18]*m(2,2);
    s.d[19] = d[46]*m(0,2) + d[64]*m(1,2) + d[19]*m(2,2);
    s.d[20] = d[47]*m(0,2) + d[65]*m(1,2) + d[20]*m(2,2);
    s.d[21] = d[48]*m(0,2) + d[66]*m(1,2) + d[21]*m(2,2);
    s.d[22] = d[49]*m(0,2) + d[67]*m(1,2) + d[22]*m(2,2);
    s.d[23] = d[50]*m(0,2) + d[68]*m(1,2) + d[23]*m(2,2);
    s.d[24] = d[51]*m(0,2) + d[69]*m(1,2) + d[24]*m(2,2);
    s.d[25] = d[52]*m(0,2) + d[70]*m(1,2) + d[25]*m(2,2);
    s.d[26] = d[53]*m(0,2) + d[71]*m(1,2) + d[26]*m(2,2);

    s.d[27] = d[0]*m(0,1) + d[27]*m(1,1) + d[72]*m(2,1);
    s.d[28] = d[1]*m(0,1) + d[28]*m(1,1) + d[73]*m(2,1);
    s.d[29] = d[2]*m(0,1) + d[29]*m(1,1) + d[74]*m(2,1);
    s.d[30] = d[3]*m(0,1) + d[30]*m(1,1) + d[75]*m(2,1);
    s.d[31] = d[4]*m(0,1) + d[31]*m(1,1) + d[76]*m(2,1);
    s.d[32] = d[5]*m(0,1) + d[32]*m(1,1) + d[77]*m(2,1);
    s.d[33] = d[6]*m(0,1) + d[33]*m(1,1) + d[78]*m(2,1);
    s.d[34] = d[7]*m(0,1) + d[34]*m(1,1) + d[79]*m(2,1);
    s.d[35] = d[8]*m(0,1) + d[35]*m(1,1) + d[80]*m(2,1);

    s.d[36] = d[54]*m(0,2) + d[ 9]*m(1,2) + d[36]*m(2,2);
    s.d[37] = d[55]*m(0,2) + d[10]*m(1,2) + d[37]*m(2,2);
    s.d[38] = d[56]*m(0,2) + d[11]*m(1,2) + d[38]*m(2,2);
    s.d[39] = d[57]*m(0,2) + d[12]*m(1,2) + d[39]*m(2,2);
    s.d[40] = d[58]*m(0,2) + d[13]*m(1,2) + d[40]*m(2,2);
    s.d[41] = d[59]*m(0,2) + d[14]*m(1,2) + d[41]*m(2,2);
    s.d[42] = d[60]*m(0,2) + d[15]*m(1,2) + d[42]*m(2,2);
    s.d[43] = d[61]*m(0,2) + d[16]*m(1,2) + d[43]*m(2,2);
    s.d[44] = d[62]*m(0,2) + d[17]*m(1,2) + d[44]*m(2,2);

    s.d[45] = d[45]*m(0,0) + d[63]*m(1,0) + d[18]*m(2,0);
    s.d[46] = d[46]*m(0,0) + d[64]*m(1,0) + d[19]*m(2,0);
    s.d[47] = d[47]*m(0,0) + d[65]*m(1,0) + d[20]*m(2,0);
    s.d[48] = d[48]*m(0,0) + d[66]*m(1,0) + d[21]*m(2,0);
    s.d[49] = d[49]*m(0,0) + d[67]*m(1,0) + d[22]*m(2,0);
    s.d[50] = d[50]*m(0,0) + d[68]*m(1,0) + d[23]*m(2,0);
    s.d[51] = d[51]*m(0,0) + d[69]*m(1,0) + d[24]*m(2,0);
    s.d[52] = d[52]*m(0,0) + d[70]*m(1,0) + d[25]*m(2,0);
    s.d[53] = d[53]*m(0,0) + d[71]*m(1,0) + d[26]*m(2,0);

    s.d[54] = d[54]*m(0,0) + d[ 9]*m(1,0) + d[36]*m(2,0);
    s.d[55] = d[55]*m(0,0) + d[10]*m(1,0) + d[37]*m(2,0);
    s.d[56] = d[56]*m(0,0) + d[11]*m(1,0) + d[38]*m(2,0);
    s.d[57] = d[57]*m(0,0) + d[12]*m(1,0) + d[39]*m(2,0);
    s.d[58] = d[58]*m(0,0) + d[13]*m(1,0) + d[40]*m(2,0);
    s.d[59] = d[59]*m(0,0) + d[14]*m(1,0) + d[41]*m(2,0);
    s.d[60] = d[60]*m(0,0) + d[15]*m(1,0) + d[42]*m(2,0);
    s.d[61] = d[61]*m(0,0) + d[16]*m(1,0) + d[43]*m(2,0);
    s.d[62] = d[62]*m(0,0) + d[17]*m(1,0) + d[44]*m(2,0);

    s.d[63] = d[45]*m(0,1) + d[63]*m(1,1) + d[18]*m(2,1);
    s.d[64] = d[46]*m(0,1) + d[64]*m(1,1) + d[19]*m(2,1);
    s.d[65] = d[47]*m(0,1) + d[65]*m(1,1) + d[20]*m(2,1);
    s.d[66] = d[48]*m(0,1) + d[66]*m(1,1) + d[21]*m(2,1);
    s.d[67] = d[49]*m(0,1) + d[67]*m(1,1) + d[22]*m(2,1);
    s.d[68] = d[50]*m(0,1) + d[68]*m(1,1) + d[23]*m(2,1);
    s.d[69] = d[51]*m(0,1) + d[69]*m(1,1) + d[24]*m(2,1);
    s.d[70] = d[52]*m(0,1) + d[70]*m(1,1) + d[25]*m(2,1);
    s.d[71] = d[53]*m(0,1) + d[71]*m(1,1) + d[26]*m(2,1);

    s.d[72] = d[0]*m(0,2) + d[27]*m(1,2) + d[72]*m(2,2);
    s.d[73] = d[1]*m(0,2) + d[28]*m(1,2) + d[73]*m(2,2);
    s.d[74] = d[2]*m(0,2) + d[29]*m(1,2) + d[74]*m(2,2);
    s.d[75] = d[3]*m(0,2) + d[30]*m(1,2) + d[75]*m(2,2);
    s.d[76] = d[4]*m(0,2) + d[31]*m(1,2) + d[76]*m(2,2);
    s.d[77] = d[5]*m(0,2) + d[32]*m(1,2) + d[77]*m(2,2);
    s.d[78] = d[6]*m(0,2) + d[33]*m(1,2) + d[78]*m(2,2);
    s.d[79] = d[7]*m(0,2) + d[34]*m(1,2) + d[79]*m(2,2);
    s.d[80] = d[8]*m(0,2) + d[35]*m(1,2) + d[80]*m(2,2);

    return s;
}

//-----------------------------------------------------------------------------
// trace
// C.tr() = I:C:I
inline double tens4d::tr() const
{
    return (d[0]+d[1]+d[2]+d[9]+d[10]+d[11]+d[18]+d[19]+d[20]);
}

//-----------------------------------------------------------------------------
// intialize to zero
inline void tens4d::zero()
{
    d[ 0] = d[ 1] = d[ 2] = d[ 3] = d[ 4] = d[ 5] = d[ 6] = d[ 7] = d[ 8] = 0;
    d[ 9] = d[10] = d[11] = d[12] = d[13] = d[14] = d[15] = d[16] = d[17] = 0;
    d[18] = d[19] = d[20] = d[21] = d[22] = d[23] = d[24] = d[25] = d[26] = 0;
    d[27] = d[28] = d[29] = d[30] = d[31] = d[32] = d[33] = d[34] = d[35] = 0;
    d[36] = d[37] = d[38] = d[39] = d[40] = d[41] = d[42] = d[43] = d[44] = 0;
    d[45] = d[46] = d[47] = d[48] = d[49] = d[50] = d[51] = d[52] = d[53] = 0;
    d[54] = d[55] = d[56] = d[57] = d[58] = d[59] = d[60] = d[61] = d[62] = 0;
    d[63] = d[64] = d[65] = d[66] = d[67] = d[68] = d[69] = d[70] = d[71] = 0;
    d[72] = d[73] = d[74] = d[75] = d[76] = d[77] = d[78] = d[79] = d[80] = 0;
}

//-----------------------------------------------------------------------------
// extract 9x9 matrix
inline void tens4d::extract(double D[9][9])
{
    D[0][0] = d[ 0];    D[0][3] = d[27];    D[0][6] = d[54];
    D[1][0] = d[ 1];    D[1][3] = d[28];    D[1][6] = d[55];
    D[2][0] = d[ 2];    D[2][3] = d[29];    D[2][6] = d[56];
    D[3][0] = d[ 3];    D[3][3] = d[30];    D[3][6] = d[57];
    D[4][0] = d[ 4];    D[4][3] = d[31];    D[4][6] = d[58];
    D[5][0] = d[ 5];    D[5][3] = d[32];    D[5][6] = d[59];
    D[6][0] = d[ 6];    D[6][3] = d[33];    D[6][6] = d[60];
    D[7][0] = d[ 7];    D[7][3] = d[34];    D[7][6] = d[61];
    D[8][0] = d[ 8];    D[8][3] = d[35];    D[8][6] = d[62];
    D[0][1] = d[ 9];    D[0][4] = d[36];    D[0][7] = d[63];
    D[1][1] = d[10];    D[1][4] = d[37];    D[1][7] = d[64];
    D[2][1] = d[11];    D[2][4] = d[38];    D[2][7] = d[65];
    D[3][1] = d[12];    D[3][4] = d[39];    D[3][7] = d[66];
    D[4][1] = d[13];    D[4][4] = d[40];    D[4][7] = d[67];
    D[5][1] = d[14];    D[5][4] = d[41];    D[5][7] = d[68];
    D[6][1] = d[15];    D[6][4] = d[42];    D[6][7] = d[69];
    D[7][1] = d[16];    D[7][4] = d[43];    D[7][7] = d[70];
    D[8][1] = d[17];    D[8][4] = d[44];    D[8][7] = d[71];
    D[0][2] = d[18];    D[0][5] = d[45];    D[0][8] = d[72];
    D[1][2] = d[19];    D[1][5] = d[46];    D[1][8] = d[73];
    D[2][2] = d[20];    D[2][5] = d[47];    D[2][8] = d[74];
    D[3][2] = d[21];    D[3][5] = d[48];    D[3][8] = d[75];
    D[4][2] = d[22];    D[4][5] = d[49];    D[4][8] = d[76];
    D[5][2] = d[23];    D[5][5] = d[50];    D[5][8] = d[77];
    D[6][2] = d[24];    D[6][5] = d[51];    D[6][8] = d[78];
    D[7][2] = d[25];    D[7][5] = d[52];    D[7][8] = d[79];
    D[8][2] = d[26];    D[8][5] = d[53];    D[8][8] = d[80];
}


//-----------------------------------------------------------------------------
// compute the super symmetric (major and minor symmetric) component of the tensor
// Sijkl = (1/8)*(Cijkl + Cjikl + Cijlk + Cjilk + Cklij + Clkij + Cklji + Clkji)
inline tens4ds tens4d::supersymm() const
{
    tens4ds s;
    
    s.d[ 0] = d[0];
    s.d[ 1] = (d[9]+d[1])/2;
    s.d[ 2] = d[10];
    s.d[ 3] = (d[18]+d[2])/2;
    s.d[ 4] = (d[19]+d[11])/2;
    s.d[ 5] = d[20];
    s.d[ 6] = (d[3] + d[6] + d[27] + d[54])/4;
    s.d[ 7] = (d[12] + d[15] + d[28] + d[55])/4;
    s.d[ 8] = (d[21] + d[24] + d[29] + d[56])/4;
    s.d[ 9] = (d[30] + d[33] + d[57] + d[60])/4;
    s.d[10] = (d[4] + d[7] + d[36] + d[63])/4;
    s.d[11] = (d[13] + d[16] + d[37] + d[64])/4;
    s.d[12] = (d[22] + d[25] + d[38] + d[65])/4;
    s.d[13] = (d[31] + d[34] + d[39] + d[42] + d[58] + d[61] + d[66] + d[69])/8;
    s.d[14] = (d[40] + d[43] + d[67] + d[70])/4;
    s.d[15] = (d[5] + d[8] + d[45] + d[72])/4;
    s.d[16] = (d[14] + d[17] + d[46] + d[73])/4;
    s.d[17] = (d[23] + d[26] + d[47] + d[74])/4;
    s.d[18] = (d[32] + d[35] + d[48] + d[51] + d[59] + d[62] + d[75] + d[78])/8;
    s.d[19] = (d[41] + d[44] + d[49] + d[52] + d[68] + d[71] + d[76] + d[79])/8;
    s.d[20] = (d[50] + d[53] + d[77] + d[80])/4;

    return s;
}

//-----------------------------------------------------------------------------
// compute the major transpose of the tensor
// Sijkl -> Sklij
inline tens4d tens4d::transpose() const
{
    tens4d s;
    
    s.d[ 0] = d[ 0];    s.d[ 6] = d[54];    s.d[ 5] = d[45];
    s.d[27] = d[ 3];    s.d[33] = d[57];    s.d[32] = d[48];
    s.d[72] = d[ 8];    s.d[78] = d[62];    s.d[77] = d[53];
    s.d[54] = d[ 6];    s.d[60] = d[60];    s.d[59] = d[51];
    s.d[ 9] = d[ 1];    s.d[15] = d[55];    s.d[14] = d[46];
    s.d[36] = d[ 4];    s.d[42] = d[58];    s.d[41] = d[49];
    s.d[45] = d[ 5];    s.d[51] = d[59];    s.d[50] = d[50];
    s.d[63] = d[ 7];    s.d[69] = d[61];    s.d[68] = d[52];
    s.d[18] = d[ 2];    s.d[24] = d[56];    s.d[23] = d[47];

    s.d[ 3] = d[27];    s.d[ 1] = d[ 9];    s.d[ 7] = d[63];
    s.d[30] = d[30];    s.d[28] = d[12];    s.d[34] = d[66];
    s.d[75] = d[35];    s.d[73] = d[17];    s.d[79] = d[71];
    s.d[57] = d[33];    s.d[55] = d[15];    s.d[61] = d[69];
    s.d[12] = d[28];    s.d[10] = d[10];    s.d[16] = d[64];
    s.d[39] = d[31];    s.d[37] = d[13];    s.d[43] = d[67];
    s.d[48] = d[32];    s.d[46] = d[14];    s.d[52] = d[68];
    s.d[66] = d[34];    s.d[64] = d[16];    s.d[70] = d[70];
    s.d[21] = d[29];    s.d[19] = d[11];    s.d[25] = d[65];

    s.d[ 8] = d[72];    s.d[ 4] = d[36];    s.d[ 2] = d[18];
    s.d[35] = d[75];    s.d[31] = d[39];    s.d[29] = d[21];
    s.d[80] = d[80];    s.d[76] = d[44];    s.d[74] = d[26];
    s.d[62] = d[78];    s.d[58] = d[42];    s.d[56] = d[24];
    s.d[17] = d[73];    s.d[13] = d[37];    s.d[11] = d[19];
    s.d[44] = d[76];    s.d[40] = d[40];    s.d[38] = d[22];
    s.d[53] = d[77];    s.d[49] = d[41];    s.d[47] = d[23];
    s.d[71] = d[79];    s.d[67] = d[43];    s.d[65] = d[25];
    s.d[26] = d[74];    s.d[22] = d[38];    s.d[20] = d[20];
    
    return s;
}

//-----------------------------------------------------------------------------
// compute the left transpose of the tensor
// Sijkl -> Sjikl
inline tens4d tens4d::left_transpose() const
{
    tens4d s;
    
    s.d[ 0] = d[ 0];    s.d[ 6] = d[ 3];    s.d[ 5] = d[ 8];
    s.d[27] = d[27];    s.d[33] = d[30];    s.d[32] = d[35];
    s.d[72] = d[72];    s.d[78] = d[75];    s.d[77] = d[80];
    s.d[54] = d[54];    s.d[60] = d[57];    s.d[59] = d[62];
    s.d[ 9] = d[ 9];    s.d[15] = d[12];    s.d[14] = d[17];
    s.d[36] = d[36];    s.d[42] = d[39];    s.d[41] = d[44];
    s.d[45] = d[45];    s.d[51] = d[48];    s.d[50] = d[53];
    s.d[63] = d[63];    s.d[69] = d[66];    s.d[68] = d[71];
    s.d[18] = d[18];    s.d[24] = d[21];    s.d[23] = d[26];

    s.d[ 3] = d[ 6];    s.d[ 1] = d[ 1];    s.d[ 7] = d[ 4];
    s.d[30] = d[33];    s.d[28] = d[28];    s.d[34] = d[31];
    s.d[75] = d[78];    s.d[73] = d[73];    s.d[79] = d[76];
    s.d[57] = d[60];    s.d[55] = d[55];    s.d[61] = d[58];
    s.d[12] = d[15];    s.d[10] = d[10];    s.d[16] = d[13];
    s.d[39] = d[42];    s.d[37] = d[37];    s.d[43] = d[40];
    s.d[48] = d[51];    s.d[46] = d[46];    s.d[52] = d[49];
    s.d[66] = d[69];    s.d[64] = d[64];    s.d[70] = d[67];
    s.d[21] = d[24];    s.d[19] = d[19];    s.d[25] = d[22];

    s.d[ 8] = d[ 5];    s.d[ 4] = d[ 7];    s.d[ 2] = d[ 2];
    s.d[35] = d[32];    s.d[31] = d[34];    s.d[29] = d[29];
    s.d[80] = d[77];    s.d[76] = d[79];    s.d[74] = d[74];
    s.d[62] = d[59];    s.d[58] = d[61];    s.d[56] = d[56];
    s.d[17] = d[14];    s.d[13] = d[16];    s.d[11] = d[11];
    s.d[44] = d[41];    s.d[40] = d[43];    s.d[38] = d[38];
    s.d[53] = d[50];    s.d[49] = d[52];    s.d[47] = d[47];
    s.d[71] = d[68];    s.d[67] = d[70];    s.d[65] = d[65];
    s.d[26] = d[23];    s.d[22] = d[25];    s.d[20] = d[20];
    
    return s;
}

//-----------------------------------------------------------------------------
// compute the right transpose of the tensor
// Sijkl -> Sijlk
inline tens4d tens4d::right_transpose() const
{
    tens4d s;
    
    s.d[ 0] = d[ 0];    s.d[ 6] = d[ 6];    s.d[ 5] = d[ 5];
    s.d[27] = d[54];    s.d[33] = d[60];    s.d[32] = d[59];
    s.d[72] = d[45];    s.d[78] = d[51];    s.d[77] = d[50];
    s.d[54] = d[27];    s.d[60] = d[33];    s.d[59] = d[32];
    s.d[ 9] = d[ 9];    s.d[15] = d[15];    s.d[14] = d[14];
    s.d[36] = d[63];    s.d[42] = d[69];    s.d[41] = d[68];
    s.d[45] = d[72];    s.d[51] = d[78];    s.d[50] = d[77];
    s.d[63] = d[36];    s.d[69] = d[42];    s.d[68] = d[41];
    s.d[18] = d[18];    s.d[24] = d[24];    s.d[23] = d[23];

    s.d[ 3] = d[ 3];    s.d[ 1] = d[ 1];    s.d[ 7] = d[ 7];
    s.d[30] = d[57];    s.d[28] = d[55];    s.d[34] = d[61];
    s.d[75] = d[48];    s.d[73] = d[46];    s.d[79] = d[52];
    s.d[57] = d[30];    s.d[55] = d[28];    s.d[61] = d[34];
    s.d[12] = d[12];    s.d[10] = d[10];    s.d[16] = d[16];
    s.d[39] = d[66];    s.d[37] = d[64];    s.d[43] = d[70];
    s.d[48] = d[75];    s.d[46] = d[73];    s.d[52] = d[79];
    s.d[66] = d[39];    s.d[64] = d[37];    s.d[70] = d[43];
    s.d[21] = d[21];    s.d[19] = d[19];    s.d[25] = d[25];

    s.d[ 8] = d[ 8];    s.d[ 4] = d[ 4];    s.d[ 2] = d[ 2];
    s.d[35] = d[62];    s.d[31] = d[58];    s.d[29] = d[56];
    s.d[80] = d[53];    s.d[76] = d[49];    s.d[74] = d[47];
    s.d[62] = d[35];    s.d[58] = d[31];    s.d[56] = d[29];
    s.d[17] = d[17];    s.d[13] = d[13];    s.d[11] = d[11];
    s.d[44] = d[71];    s.d[40] = d[67];    s.d[38] = d[65];
    s.d[53] = d[80];    s.d[49] = d[76];    s.d[47] = d[74];
    s.d[71] = d[44];    s.d[67] = d[40];    s.d[65] = d[38];
    s.d[26] = d[26];    s.d[22] = d[22];    s.d[20] = d[20];
    
    return s;
}

//-----------------------------------------------------------------------------
// (a dyad1 b)_ijkl = a_ij b_kl
inline tens4d dyad1(const mat3d& a, const mat3d& b)
{
    tens4d c;
    
    c.d[ 0] = a(0,0)*b(0,0);
    c.d[27] = a(0,0)*b(0,1);
    c.d[72] = a(0,0)*b(0,2);
    c.d[54] = a(0,0)*b(1,0);
    c.d[ 9] = a(0,0)*b(1,1);
    c.d[36] = a(0,0)*b(1,2);
    c.d[45] = a(0,0)*b(2,0);
    c.d[63] = a(0,0)*b(2,1);
    c.d[18] = a(0,0)*b(2,2);

    c.d[ 3] = a(0,1)*b(0,0);
    c.d[30] = a(0,1)*b(0,1);
    c.d[75] = a(0,1)*b(0,2);
    c.d[57] = a(0,1)*b(1,0);
    c.d[12] = a(0,1)*b(1,1);
    c.d[39] = a(0,1)*b(1,2);
    c.d[48] = a(0,1)*b(2,0);
    c.d[66] = a(0,1)*b(2,1);
    c.d[21] = a(0,1)*b(2,2);

    c.d[ 8] = a(0,2)*b(0,0);
    c.d[35] = a(0,2)*b(0,1);
    c.d[80] = a(0,2)*b(0,2);
    c.d[62] = a(0,2)*b(1,0);
    c.d[17] = a(0,2)*b(1,1);
    c.d[44] = a(0,2)*b(1,2);
    c.d[53] = a(0,2)*b(2,0);
    c.d[71] = a(0,2)*b(2,1);
    c.d[26] = a(0,2)*b(2,2);

    c.d[6] = a(1,0)*b(0,0);
    c.d[33] = a(1,0)*b(0,1);
    c.d[78] = a(1,0)*b(0,2);
    c.d[60] = a(1,0)*b(1,0);
    c.d[15] = a(1,0)*b(1,1);
    c.d[42] = a(1,0)*b(1,2);
    c.d[51] = a(1,0)*b(2,0);
    c.d[69] = a(1,0)*b(2,1);
    c.d[24] = a(1,0)*b(2,2);

    c.d[ 1] = a(1,1)*b(0,0);
    c.d[28] = a(1,1)*b(0,1);
    c.d[73] = a(1,1)*b(0,2);
    c.d[55] = a(1,1)*b(1,0);
    c.d[10] = a(1,1)*b(1,1);
    c.d[37] = a(1,1)*b(1,2);
    c.d[46] = a(1,1)*b(2,0);
    c.d[64] = a(1,1)*b(2,1);
    c.d[19] = a(1,1)*b(2,2);

    c.d[ 4] = a(1,2)*b(0,0);
    c.d[31] = a(1,2)*b(0,1);
    c.d[76] = a(1,2)*b(0,2);
    c.d[58] = a(1,2)*b(1,0);
    c.d[13] = a(1,2)*b(1,1);
    c.d[40] = a(1,2)*b(1,2);
    c.d[49] = a(1,2)*b(2,0);
    c.d[67] = a(1,2)*b(2,1);
    c.d[22] = a(1,2)*b(2,2);

    c.d[ 5] = a(2,0)*b(0,0);
    c.d[32] = a(2,0)*b(0,1);
    c.d[77] = a(2,0)*b(0,2);
    c.d[59] = a(2,0)*b(1,0);
    c.d[14] = a(2,0)*b(1,1);
    c.d[41] = a(2,0)*b(1,2);
    c.d[50] = a(2,0)*b(2,0);
    c.d[68] = a(2,0)*b(2,1);
    c.d[23] = a(2,0)*b(2,2);

    c.d[ 7] = a(2,1)*b(0,0);
    c.d[34] = a(2,1)*b(0,1);
    c.d[79] = a(2,1)*b(0,2);
    c.d[61] = a(2,1)*b(1,0);
    c.d[16] = a(2,1)*b(1,1);
    c.d[43] = a(2,1)*b(1,2);
    c.d[52] = a(2,1)*b(2,0);
    c.d[70] = a(2,1)*b(2,1);
    c.d[25] = a(2,1)*b(2,2);

    c.d[ 2] = a(2,2)*b(0,0);
    c.d[29] = a(2,2)*b(0,1);
    c.d[74] = a(2,2)*b(0,2);
    c.d[56] = a(2,2)*b(1,0);
    c.d[11] = a(2,2)*b(1,1);
    c.d[38] = a(2,2)*b(1,2);
    c.d[47] = a(2,2)*b(2,0);
    c.d[65] = a(2,2)*b(2,1);
    c.d[20] = a(2,2)*b(2,2);

    return c;
}

//-----------------------------------------------------------------------------
// (a dyad2 b)_ijkl = a_ik b_jl
inline tens4d dyad2(const mat3d& a, const mat3d& b)
{
    tens4d c;
    
    c.d[ 0] = a(0,0)*b(0,0);
    c.d[27] = a(0,0)*b(0,1);
    c.d[72] = a(0,0)*b(0,2);
    c.d[54] = a(0,1)*b(0,0);
    c.d[ 9] = a(0,1)*b(0,1);
    c.d[36] = a(0,1)*b(0,2);
    c.d[45] = a(0,2)*b(0,0);
    c.d[63] = a(0,2)*b(0,1);
    c.d[18] = a(0,2)*b(0,2);

    c.d[ 3] = a(0,0)*b(1,0);
    c.d[30] = a(0,0)*b(1,1);
    c.d[75] = a(0,0)*b(1,2);
    c.d[57] = a(0,1)*b(1,0);
    c.d[12] = a(0,1)*b(1,1);
    c.d[39] = a(0,1)*b(1,2);
    c.d[48] = a(0,2)*b(1,0);
    c.d[66] = a(0,2)*b(1,1);
    c.d[21] = a(0,2)*b(1,2);

    c.d[ 8] = a(0,0)*b(2,0);
    c.d[35] = a(0,0)*b(2,1);
    c.d[80] = a(0,0)*b(2,2);
    c.d[62] = a(0,1)*b(2,0);
    c.d[17] = a(0,1)*b(2,1);
    c.d[44] = a(0,1)*b(2,2);
    c.d[53] = a(0,2)*b(2,0);
    c.d[71] = a(0,2)*b(2,1);
    c.d[26] = a(0,2)*b(2,2);

    c.d[ 6] = a(1,0)*b(0,0);
    c.d[33] = a(1,0)*b(0,1);
    c.d[78] = a(1,0)*b(0,2);
    c.d[60] = a(1,1)*b(0,0);
    c.d[15] = a(1,1)*b(0,1);
    c.d[42] = a(1,1)*b(0,2);
    c.d[51] = a(1,2)*b(0,0);
    c.d[69] = a(1,2)*b(0,1);
    c.d[24] = a(1,2)*b(0,2);

    c.d[ 1] = a(1,0)*b(1,0);
    c.d[28] = a(1,0)*b(1,1);
    c.d[73] = a(1,0)*b(1,2);
    c.d[55] = a(1,1)*b(1,0);
    c.d[10] = a(1,1)*b(1,1);
    c.d[37] = a(1,1)*b(1,2);
    c.d[46] = a(1,2)*b(1,0);
    c.d[64] = a(1,2)*b(1,1);
    c.d[19] = a(1,2)*b(1,2);

    c.d[ 4] = a(1,0)*b(2,0);
    c.d[31] = a(1,0)*b(2,1);
    c.d[76] = a(1,0)*b(2,2);
    c.d[58] = a(1,1)*b(2,0);
    c.d[13] = a(1,1)*b(2,1);
    c.d[40] = a(1,1)*b(2,2);
    c.d[49] = a(1,2)*b(2,0);
    c.d[67] = a(1,2)*b(2,1);
    c.d[22] = a(1,2)*b(2,2);

    c.d[ 5] = a(2,0)*b(0,0);
    c.d[32] = a(2,0)*b(0,1);
    c.d[77] = a(2,0)*b(0,2);
    c.d[59] = a(2,1)*b(0,0);
    c.d[14] = a(2,1)*b(0,1);
    c.d[41] = a(2,1)*b(0,2);
    c.d[50] = a(2,2)*b(0,0);
    c.d[68] = a(2,2)*b(0,1);
    c.d[23] = a(2,2)*b(0,2);

    c.d[ 7] = a(2,0)*b(1,0);
    c.d[34] = a(2,0)*b(1,1);
    c.d[79] = a(2,0)*b(1,2);
    c.d[61] = a(2,1)*b(1,0);
    c.d[16] = a(2,1)*b(1,1);
    c.d[43] = a(2,1)*b(1,2);
    c.d[52] = a(2,2)*b(1,0);
    c.d[70] = a(2,2)*b(1,1);
    c.d[25] = a(2,2)*b(1,2);

    c.d[ 2] = a(2,0)*b(2,0);
    c.d[29] = a(2,0)*b(2,1);
    c.d[74] = a(2,0)*b(2,2);
    c.d[56] = a(2,1)*b(2,0);
    c.d[11] = a(2,1)*b(2,1);
    c.d[38] = a(2,1)*b(2,2);
    c.d[47] = a(2,2)*b(2,0);
    c.d[65] = a(2,2)*b(2,1);
    c.d[20] = a(2,2)*b(2,2);

    return c;
}

//-----------------------------------------------------------------------------
// (a dyad3 b)_ijkl = a_il b_jk
inline tens4d dyad3(const mat3d& a, const mat3d& b)
{
    tens4d c;
    
    c.d[ 0] = a(0,0)*b(0,0);
    c.d[27] = a(0,1)*b(0,0);
    c.d[72] = a(0,2)*b(0,0);
    c.d[54] = a(0,0)*b(0,1);
    c.d[ 9] = a(0,1)*b(0,1);
    c.d[36] = a(0,2)*b(0,1);
    c.d[45] = a(0,0)*b(0,2);
    c.d[63] = a(0,1)*b(0,2);
    c.d[18] = a(0,2)*b(0,2);

    c.d[ 3] = a(0,0)*b(1,0);
    c.d[30] = a(0,1)*b(1,0);
    c.d[75] = a(0,2)*b(1,0);
    c.d[57] = a(0,0)*b(1,1);
    c.d[12] = a(0,1)*b(1,1);
    c.d[39] = a(0,2)*b(1,1);
    c.d[48] = a(0,0)*b(1,2);
    c.d[66] = a(0,1)*b(1,2);
    c.d[21] = a(0,2)*b(1,2);

    c.d[ 8] = a(0,0)*b(2,0);
    c.d[35] = a(0,1)*b(2,0);
    c.d[80] = a(0,2)*b(2,0);
    c.d[62] = a(0,0)*b(2,1);
    c.d[17] = a(0,1)*b(2,1);
    c.d[44] = a(0,2)*b(2,1);
    c.d[53] = a(0,0)*b(2,2);
    c.d[71] = a(0,1)*b(2,2);
    c.d[26] = a(0,2)*b(2,2);

    c.d[ 6] = a(1,0)*b(0,0);
    c.d[33] = a(1,1)*b(0,0);
    c.d[78] = a(1,2)*b(0,0);
    c.d[60] = a(1,0)*b(0,1);
    c.d[15] = a(1,1)*b(0,1);
    c.d[42] = a(1,2)*b(0,1);
    c.d[51] = a(1,0)*b(0,2);
    c.d[69] = a(1,1)*b(0,2);
    c.d[24] = a(1,2)*b(0,2);

    c.d[ 1] = a(1,0)*b(1,0);
    c.d[28] = a(1,1)*b(1,0);
    c.d[73] = a(1,2)*b(1,0);
    c.d[55] = a(1,0)*b(1,1);
    c.d[10] = a(1,1)*b(1,1);
    c.d[37] = a(1,2)*b(1,1);
    c.d[46] = a(1,0)*b(1,2);
    c.d[64] = a(1,1)*b(1,2);
    c.d[19] = a(1,2)*b(1,2);

    c.d[ 4] = a(1,0)*b(2,0);
    c.d[31] = a(1,1)*b(2,0);
    c.d[76] = a(1,2)*b(2,0);
    c.d[58] = a(1,0)*b(2,1);
    c.d[13] = a(1,1)*b(2,1);
    c.d[40] = a(1,2)*b(2,1);
    c.d[49] = a(1,0)*b(2,2);
    c.d[67] = a(1,1)*b(2,2);
    c.d[22] = a(1,2)*b(2,2);

    c.d[ 5] = a(2,0)*b(0,0);
    c.d[32] = a(2,1)*b(0,0);
    c.d[77] = a(2,2)*b(0,0);
    c.d[59] = a(2,0)*b(0,1);
    c.d[14] = a(2,1)*b(0,1);
    c.d[41] = a(2,2)*b(0,1);
    c.d[50] = a(2,0)*b(0,2);
    c.d[68] = a(2,1)*b(0,2);
    c.d[23] = a(2,2)*b(0,2);

    c.d[ 7] = a(2,0)*b(1,0);
    c.d[34] = a(2,1)*b(1,0);
    c.d[79] = a(2,2)*b(1,0);
    c.d[61] = a(2,0)*b(1,1);
    c.d[16] = a(2,1)*b(1,1);
    c.d[43] = a(2,2)*b(1,1);
    c.d[52] = a(2,0)*b(1,2);
    c.d[70] = a(2,1)*b(1,2);
    c.d[25] = a(2,2)*b(1,2);

    c.d[ 2] = a(2,0)*b(2,0);
    c.d[29] = a(2,1)*b(2,0);
    c.d[74] = a(2,2)*b(2,0);
    c.d[56] = a(2,0)*b(2,1);
    c.d[11] = a(2,1)*b(2,1);
    c.d[38] = a(2,2)*b(2,1);
    c.d[47] = a(2,0)*b(2,2);
    c.d[65] = a(2,1)*b(2,2);
    c.d[20] = a(2,2)*b(2,2);

    return c;
}

//-----------------------------------------------------------------------------
// (a ddot b)_ijkl = a_ijmn b_mnkl
inline tens4d ddot(const tens4d& a, const tens4d& b)
{
    tens4d c;
    
    c.d[0] = a.d[0]*b.d[0] + a.d[9]*b.d[1] + a.d[18]*b.d[2] + a.d[27]*b.d[3] + a.d[36]*b.d[4] + a.d[45]*b.d[5] + a.d[54]*b.d[6] + a.d[63]*b.d[7] + a.d[72]*b.d[8];
    c.d[1] = a.d[0]*b.d[9] + a.d[9]*b.d[10] + a.d[18]*b.d[11] + a.d[27]*b.d[12] + a.d[36]*b.d[13] + a.d[45]*b.d[14] + a.d[54]*b.d[15] + a.d[63]*b.d[16] + a.d[72]*b.d[17];
    c.d[2] = a.d[0]*b.d[18] + a.d[9]*b.d[19] + a.d[18]*b.d[20] + a.d[27]*b.d[21] + a.d[36]*b.d[22] + a.d[45]*b.d[23] + a.d[54]*b.d[24] + a.d[63]*b.d[25] + a.d[72]*b.d[26];
    c.d[3] = a.d[0]*b.d[27] + a.d[9]*b.d[28] + a.d[18]*b.d[29] + a.d[27]*b.d[30] + a.d[36]*b.d[31] + a.d[45]*b.d[32] + a.d[54]*b.d[33] + a.d[63]*b.d[34] + a.d[72]*b.d[35];
    c.d[4] = a.d[0]*b.d[36] + a.d[9]*b.d[37] + a.d[18]*b.d[38] + a.d[27]*b.d[39] + a.d[36]*b.d[40] + a.d[45]*b.d[41] + a.d[54]*b.d[42] + a.d[63]*b.d[43] + a.d[72]*b.d[44];
    c.d[5] = a.d[0]*b.d[45] + a.d[9]*b.d[46] + a.d[18]*b.d[47] + a.d[27]*b.d[48] + a.d[36]*b.d[49] + a.d[45]*b.d[50] + a.d[54]*b.d[51] + a.d[63]*b.d[52] + a.d[72]*b.d[53];
    c.d[6] = a.d[0]*b.d[54] + a.d[9]*b.d[55] + a.d[18]*b.d[56] + a.d[27]*b.d[57] + a.d[36]*b.d[58] + a.d[45]*b.d[59] + a.d[54]*b.d[60] + a.d[63]*b.d[61] + a.d[72]*b.d[62];
    c.d[7] = a.d[0]*b.d[63] + a.d[9]*b.d[64] + a.d[18]*b.d[65] + a.d[27]*b.d[66] + a.d[36]*b.d[67] + a.d[45]*b.d[68] + a.d[54]*b.d[69] + a.d[63]*b.d[70] + a.d[72]*b.d[71];
    c.d[8] = a.d[0]*b.d[72] + a.d[9]*b.d[73] + a.d[18]*b.d[74] + a.d[27]*b.d[75] + a.d[36]*b.d[76] + a.d[45]*b.d[77] + a.d[54]*b.d[78] + a.d[63]*b.d[79] + a.d[72]*b.d[80];

    c.d[9] = a.d[1]*b.d[0] + a.d[10]*b.d[1] + a.d[19]*b.d[2] + a.d[28]*b.d[3] + a.d[37]*b.d[4] + a.d[46]*b.d[5] + a.d[55]*b.d[6] + a.d[64]*b.d[7] + a.d[73]*b.d[8];
    c.d[10] = a.d[1]*b.d[9] + a.d[10]*b.d[10] + a.d[19]*b.d[11] + a.d[28]*b.d[12] + a.d[37]*b.d[13] + a.d[46]*b.d[14] + a.d[55]*b.d[15] + a.d[64]*b.d[16] + a.d[73]*b.d[17];
    c.d[11] = a.d[1]*b.d[18] + a.d[10]*b.d[19] + a.d[19]*b.d[20] + a.d[28]*b.d[21] + a.d[37]*b.d[22] + a.d[46]*b.d[23] + a.d[55]*b.d[24] + a.d[64]*b.d[25] + a.d[73]*b.d[26];
    c.d[12] = a.d[1]*b.d[27] + a.d[10]*b.d[28] + a.d[19]*b.d[29] + a.d[28]*b.d[30] + a.d[37]*b.d[31] + a.d[46]*b.d[32] + a.d[55]*b.d[33] + a.d[64]*b.d[34] + a.d[73]*b.d[35];
    c.d[13] = a.d[1]*b.d[36] + a.d[10]*b.d[37] + a.d[19]*b.d[38] + a.d[28]*b.d[39] + a.d[37]*b.d[40] + a.d[46]*b.d[41] + a.d[55]*b.d[42] + a.d[64]*b.d[43] + a.d[73]*b.d[44];
    c.d[14] = a.d[1]*b.d[45] + a.d[10]*b.d[46] + a.d[19]*b.d[47] + a.d[28]*b.d[48] + a.d[37]*b.d[49] + a.d[46]*b.d[50] + a.d[55]*b.d[51] + a.d[64]*b.d[52] + a.d[73]*b.d[53];
    c.d[15] = a.d[1]*b.d[54] + a.d[10]*b.d[55] + a.d[19]*b.d[56] + a.d[28]*b.d[57] + a.d[37]*b.d[58] + a.d[46]*b.d[59] + a.d[55]*b.d[60] + a.d[64]*b.d[61] + a.d[73]*b.d[62];
    c.d[16] = a.d[1]*b.d[63] + a.d[10]*b.d[64] + a.d[19]*b.d[65] + a.d[28]*b.d[66] + a.d[37]*b.d[67] + a.d[46]*b.d[68] + a.d[55]*b.d[69] + a.d[64]*b.d[70] + a.d[73]*b.d[71];
    c.d[17] = a.d[1]*b.d[72] + a.d[10]*b.d[73] + a.d[19]*b.d[74] + a.d[28]*b.d[75] + a.d[37]*b.d[76] + a.d[46]*b.d[77] + a.d[55]*b.d[78] + a.d[64]*b.d[79] + a.d[73]*b.d[80];

    c.d[18] = a.d[2]*b.d[0] + a.d[11]*b.d[1] + a.d[20]*b.d[2] + a.d[29]*b.d[3] + a.d[38]*b.d[4] + a.d[47]*b.d[5] + a.d[56]*b.d[6] + a.d[65]*b.d[7] + a.d[74]*b.d[8];
    c.d[19] = a.d[2]*b.d[9] + a.d[11]*b.d[10] + a.d[20]*b.d[11] + a.d[29]*b.d[12] + a.d[38]*b.d[13] + a.d[47]*b.d[14] + a.d[56]*b.d[15] + a.d[65]*b.d[16] + a.d[74]*b.d[17];
    c.d[20] = a.d[2]*b.d[18] + a.d[11]*b.d[19] + a.d[20]*b.d[20] + a.d[29]*b.d[21] + a.d[38]*b.d[22] + a.d[47]*b.d[23] + a.d[56]*b.d[24] + a.d[65]*b.d[25] + a.d[74]*b.d[26];
    c.d[21] = a.d[2]*b.d[27] + a.d[11]*b.d[28] + a.d[20]*b.d[29] + a.d[29]*b.d[30] + a.d[38]*b.d[31] + a.d[47]*b.d[32] + a.d[56]*b.d[33] + a.d[65]*b.d[34] + a.d[74]*b.d[35];
    c.d[22] = a.d[2]*b.d[36] + a.d[11]*b.d[37] + a.d[20]*b.d[38] + a.d[29]*b.d[39] + a.d[38]*b.d[40] + a.d[47]*b.d[41] + a.d[56]*b.d[42] + a.d[65]*b.d[43] + a.d[74]*b.d[44];
    c.d[23] = a.d[2]*b.d[45] + a.d[11]*b.d[46] + a.d[20]*b.d[47] + a.d[29]*b.d[48] + a.d[38]*b.d[49] + a.d[47]*b.d[50] + a.d[56]*b.d[51] + a.d[65]*b.d[52] + a.d[74]*b.d[53];
    c.d[24] = a.d[2]*b.d[54] + a.d[11]*b.d[55] + a.d[20]*b.d[56] + a.d[29]*b.d[57] + a.d[38]*b.d[58] + a.d[47]*b.d[59] + a.d[56]*b.d[60] + a.d[65]*b.d[61] + a.d[74]*b.d[62];
    c.d[25] = a.d[2]*b.d[63] + a.d[11]*b.d[64] + a.d[20]*b.d[65] + a.d[29]*b.d[66] + a.d[38]*b.d[67] + a.d[47]*b.d[68] + a.d[56]*b.d[69] + a.d[65]*b.d[70] + a.d[74]*b.d[71];
    c.d[26] = a.d[2]*b.d[72] + a.d[11]*b.d[73] + a.d[20]*b.d[74] + a.d[29]*b.d[75] + a.d[38]*b.d[76] + a.d[47]*b.d[77] + a.d[56]*b.d[78] + a.d[65]*b.d[79] + a.d[74]*b.d[80];

    c.d[27] = a.d[3]*b.d[0] + a.d[12]*b.d[1] + a.d[21]*b.d[2] + a.d[30]*b.d[3] + a.d[39]*b.d[4] + a.d[48]*b.d[5] + a.d[57]*b.d[6] + a.d[66]*b.d[7] + a.d[75]*b.d[8];
    c.d[28] = a.d[3]*b.d[9] + a.d[12]*b.d[10] + a.d[21]*b.d[11] + a.d[30]*b.d[12] + a.d[39]*b.d[13] + a.d[48]*b.d[14] + a.d[57]*b.d[15] + a.d[66]*b.d[16] + a.d[75]*b.d[17];
    c.d[29] = a.d[3]*b.d[18] + a.d[12]*b.d[19] + a.d[21]*b.d[20] + a.d[30]*b.d[21] + a.d[39]*b.d[22] + a.d[48]*b.d[23] + a.d[57]*b.d[24] + a.d[66]*b.d[25] + a.d[75]*b.d[26];
    c.d[30] = a.d[3]*b.d[27] + a.d[12]*b.d[28] + a.d[21]*b.d[29] + a.d[30]*b.d[30] + a.d[39]*b.d[31] + a.d[48]*b.d[32] + a.d[57]*b.d[33] + a.d[66]*b.d[34] + a.d[75]*b.d[35];
    c.d[31] = a.d[3]*b.d[36] + a.d[12]*b.d[37] + a.d[21]*b.d[38] + a.d[30]*b.d[39] + a.d[39]*b.d[40] + a.d[48]*b.d[41] + a.d[57]*b.d[42] + a.d[66]*b.d[43] + a.d[75]*b.d[44];
    c.d[32] = a.d[3]*b.d[45] + a.d[12]*b.d[46] + a.d[21]*b.d[47] + a.d[30]*b.d[48] + a.d[39]*b.d[49] + a.d[48]*b.d[50] + a.d[57]*b.d[51] + a.d[66]*b.d[52] + a.d[75]*b.d[53];
    c.d[33] = a.d[3]*b.d[54] + a.d[12]*b.d[55] + a.d[21]*b.d[56] + a.d[30]*b.d[57] + a.d[39]*b.d[58] + a.d[48]*b.d[59] + a.d[57]*b.d[60] + a.d[66]*b.d[61] + a.d[75]*b.d[62];
    c.d[34] = a.d[3]*b.d[63] + a.d[12]*b.d[64] + a.d[21]*b.d[65] + a.d[30]*b.d[66] + a.d[39]*b.d[67] + a.d[48]*b.d[68] + a.d[57]*b.d[69] + a.d[66]*b.d[70] + a.d[75]*b.d[71];
    c.d[35] = a.d[3]*b.d[72] + a.d[12]*b.d[73] + a.d[21]*b.d[74] + a.d[30]*b.d[75] + a.d[39]*b.d[76] + a.d[48]*b.d[77] + a.d[57]*b.d[78] + a.d[66]*b.d[79] + a.d[75]*b.d[80];

    c.d[36] = a.d[4]*b.d[0] + a.d[13]*b.d[1] + a.d[22]*b.d[2] + a.d[31]*b.d[3] + a.d[40]*b.d[4] + a.d[49]*b.d[5] + a.d[58]*b.d[6] + a.d[67]*b.d[7] + a.d[76]*b.d[8];
    c.d[37] = a.d[4]*b.d[9] + a.d[13]*b.d[10] + a.d[22]*b.d[11] + a.d[31]*b.d[12] + a.d[40]*b.d[13] + a.d[49]*b.d[14] + a.d[58]*b.d[15] + a.d[67]*b.d[16] + a.d[76]*b.d[17];
    c.d[38] = a.d[4]*b.d[18] + a.d[13]*b.d[19] + a.d[22]*b.d[20] + a.d[31]*b.d[21] + a.d[40]*b.d[22] + a.d[49]*b.d[23] + a.d[58]*b.d[24] + a.d[67]*b.d[25] + a.d[76]*b.d[26];
    c.d[39] = a.d[4]*b.d[27] + a.d[13]*b.d[28] + a.d[22]*b.d[29] + a.d[31]*b.d[30] + a.d[40]*b.d[31] + a.d[49]*b.d[32] + a.d[58]*b.d[33] + a.d[67]*b.d[34] + a.d[76]*b.d[35];
    c.d[40] = a.d[4]*b.d[36] + a.d[13]*b.d[37] + a.d[22]*b.d[38] + a.d[31]*b.d[39] + a.d[40]*b.d[40] + a.d[49]*b.d[41] + a.d[58]*b.d[42] + a.d[67]*b.d[43] + a.d[76]*b.d[44];
    c.d[41] = a.d[4]*b.d[45] + a.d[13]*b.d[46] + a.d[22]*b.d[47] + a.d[31]*b.d[48] + a.d[40]*b.d[49] + a.d[49]*b.d[50] + a.d[58]*b.d[51] + a.d[67]*b.d[52] + a.d[76]*b.d[53];
    c.d[42] = a.d[4]*b.d[54] + a.d[13]*b.d[55] + a.d[22]*b.d[56] + a.d[31]*b.d[57] + a.d[40]*b.d[58] + a.d[49]*b.d[59] + a.d[58]*b.d[60] + a.d[67]*b.d[61] + a.d[76]*b.d[62];
    c.d[43] = a.d[4]*b.d[63] + a.d[13]*b.d[64] + a.d[22]*b.d[65] + a.d[31]*b.d[66] + a.d[40]*b.d[67] + a.d[49]*b.d[68] + a.d[58]*b.d[69] + a.d[67]*b.d[70] + a.d[76]*b.d[71];
    c.d[44] = a.d[4]*b.d[72] + a.d[13]*b.d[73] + a.d[22]*b.d[74] + a.d[31]*b.d[75] + a.d[40]*b.d[76] + a.d[49]*b.d[77] + a.d[58]*b.d[78] + a.d[67]*b.d[79] + a.d[76]*b.d[80];

    c.d[45] = a.d[5]*b.d[0] + a.d[14]*b.d[1] + a.d[23]*b.d[2] + a.d[32]*b.d[3] + a.d[41]*b.d[4] + a.d[50]*b.d[5] + a.d[59]*b.d[6] + a.d[68]*b.d[7] + a.d[77]*b.d[8];
    c.d[46] = a.d[5]*b.d[9] + a.d[14]*b.d[10] + a.d[23]*b.d[11] + a.d[32]*b.d[12] + a.d[41]*b.d[13] + a.d[50]*b.d[14] + a.d[59]*b.d[15] + a.d[68]*b.d[16] + a.d[77]*b.d[17];
    c.d[47] = a.d[5]*b.d[18] + a.d[14]*b.d[19] + a.d[23]*b.d[20] + a.d[32]*b.d[21] + a.d[41]*b.d[22] + a.d[50]*b.d[23] + a.d[59]*b.d[24] + a.d[68]*b.d[25] + a.d[77]*b.d[26];
    c.d[48] = a.d[5]*b.d[27] + a.d[14]*b.d[28] + a.d[23]*b.d[29] + a.d[32]*b.d[30] + a.d[41]*b.d[31] + a.d[50]*b.d[32] + a.d[59]*b.d[33] + a.d[68]*b.d[34] + a.d[77]*b.d[35];
    c.d[49] = a.d[5]*b.d[36] + a.d[14]*b.d[37] + a.d[23]*b.d[38] + a.d[32]*b.d[39] + a.d[41]*b.d[40] + a.d[50]*b.d[41] + a.d[59]*b.d[42] + a.d[68]*b.d[43] + a.d[77]*b.d[44];
    c.d[50] = a.d[5]*b.d[45] + a.d[14]*b.d[46] + a.d[23]*b.d[47] + a.d[32]*b.d[48] + a.d[41]*b.d[49] + a.d[50]*b.d[50] + a.d[59]*b.d[51] + a.d[68]*b.d[52] + a.d[77]*b.d[53];
    c.d[51] = a.d[5]*b.d[54] + a.d[14]*b.d[55] + a.d[23]*b.d[56] + a.d[32]*b.d[57] + a.d[41]*b.d[58] + a.d[50]*b.d[59] + a.d[59]*b.d[60] + a.d[68]*b.d[61] + a.d[77]*b.d[62];
    c.d[52] = a.d[5]*b.d[63] + a.d[14]*b.d[64] + a.d[23]*b.d[65] + a.d[32]*b.d[66] + a.d[41]*b.d[67] + a.d[50]*b.d[68] + a.d[59]*b.d[69] + a.d[68]*b.d[70] + a.d[77]*b.d[71];
    c.d[53] = a.d[5]*b.d[72] + a.d[14]*b.d[73] + a.d[23]*b.d[74] + a.d[32]*b.d[75] + a.d[41]*b.d[76] + a.d[50]*b.d[77] + a.d[59]*b.d[78] + a.d[68]*b.d[79] + a.d[77]*b.d[80];

    c.d[54] = a.d[6]*b.d[0] + a.d[15]*b.d[1] + a.d[24]*b.d[2] + a.d[33]*b.d[3] + a.d[42]*b.d[4] + a.d[51]*b.d[5] + a.d[60]*b.d[6] + a.d[69]*b.d[7] + a.d[78]*b.d[8];
    c.d[55] = a.d[6]*b.d[9] + a.d[15]*b.d[10] + a.d[24]*b.d[11] + a.d[33]*b.d[12] + a.d[42]*b.d[13] + a.d[51]*b.d[14] + a.d[60]*b.d[15] + a.d[69]*b.d[16] + a.d[78]*b.d[17];
    c.d[56] = a.d[6]*b.d[18] + a.d[15]*b.d[19] + a.d[24]*b.d[20] + a.d[33]*b.d[21] + a.d[42]*b.d[22] + a.d[51]*b.d[23] + a.d[60]*b.d[24] + a.d[69]*b.d[25] + a.d[78]*b.d[26];
    c.d[57] = a.d[6]*b.d[27] + a.d[15]*b.d[28] + a.d[24]*b.d[29] + a.d[33]*b.d[30] + a.d[42]*b.d[31] + a.d[51]*b.d[32] + a.d[60]*b.d[33] + a.d[69]*b.d[34] + a.d[78]*b.d[35];
    c.d[58] = a.d[6]*b.d[36] + a.d[15]*b.d[37] + a.d[24]*b.d[38] + a.d[33]*b.d[39] + a.d[42]*b.d[40] + a.d[51]*b.d[41] + a.d[60]*b.d[42] + a.d[69]*b.d[43] + a.d[78]*b.d[44];
    c.d[59] = a.d[6]*b.d[45] + a.d[15]*b.d[46] + a.d[24]*b.d[47] + a.d[33]*b.d[48] + a.d[42]*b.d[49] + a.d[51]*b.d[50] + a.d[60]*b.d[51] + a.d[69]*b.d[52] + a.d[78]*b.d[53];
    c.d[60] = a.d[6]*b.d[54] + a.d[15]*b.d[55] + a.d[24]*b.d[56] + a.d[33]*b.d[57] + a.d[42]*b.d[58] + a.d[51]*b.d[59] + a.d[60]*b.d[60] + a.d[69]*b.d[61] + a.d[78]*b.d[62];
    c.d[61] = a.d[6]*b.d[63] + a.d[15]*b.d[64] + a.d[24]*b.d[65] + a.d[33]*b.d[66] + a.d[42]*b.d[67] + a.d[51]*b.d[68] + a.d[60]*b.d[69] + a.d[69]*b.d[70] + a.d[78]*b.d[71];
    c.d[62] = a.d[6]*b.d[72] + a.d[15]*b.d[73] + a.d[24]*b.d[74] + a.d[33]*b.d[75] + a.d[42]*b.d[76] + a.d[51]*b.d[77] + a.d[60]*b.d[78] + a.d[69]*b.d[79] + a.d[78]*b.d[80];

    c.d[63] = a.d[7]*b.d[0] + a.d[16]*b.d[1] + a.d[25]*b.d[2] + a.d[34]*b.d[3] + a.d[43]*b.d[4] + a.d[52]*b.d[5] + a.d[61]*b.d[6] + a.d[70]*b.d[7] + a.d[79]*b.d[8];
    c.d[64] = a.d[7]*b.d[9] + a.d[16]*b.d[10] + a.d[25]*b.d[11] + a.d[34]*b.d[12] + a.d[43]*b.d[13] + a.d[52]*b.d[14] + a.d[61]*b.d[15] + a.d[70]*b.d[16] + a.d[79]*b.d[17];
    c.d[65] = a.d[7]*b.d[18] + a.d[16]*b.d[19] + a.d[25]*b.d[20] + a.d[34]*b.d[21] + a.d[43]*b.d[22] + a.d[52]*b.d[23] + a.d[61]*b.d[24] + a.d[70]*b.d[25] + a.d[79]*b.d[26];
    c.d[66] = a.d[7]*b.d[27] + a.d[16]*b.d[28] + a.d[25]*b.d[29] + a.d[34]*b.d[30] + a.d[43]*b.d[31] + a.d[52]*b.d[32] + a.d[61]*b.d[33] + a.d[70]*b.d[34] + a.d[79]*b.d[35];
    c.d[67] = a.d[7]*b.d[36] + a.d[16]*b.d[37] + a.d[25]*b.d[38] + a.d[34]*b.d[39] + a.d[43]*b.d[40] + a.d[52]*b.d[41] + a.d[61]*b.d[42] + a.d[70]*b.d[43] + a.d[79]*b.d[44];
    c.d[68] = a.d[7]*b.d[45] + a.d[16]*b.d[46] + a.d[25]*b.d[47] + a.d[34]*b.d[48] + a.d[43]*b.d[49] + a.d[52]*b.d[50] + a.d[61]*b.d[51] + a.d[70]*b.d[52] + a.d[79]*b.d[53];
    c.d[69] = a.d[7]*b.d[54] + a.d[16]*b.d[55] + a.d[25]*b.d[56] + a.d[34]*b.d[57] + a.d[43]*b.d[58] + a.d[52]*b.d[59] + a.d[61]*b.d[60] + a.d[70]*b.d[61] + a.d[79]*b.d[62];
    c.d[70] = a.d[7]*b.d[63] + a.d[16]*b.d[64] + a.d[25]*b.d[65] + a.d[34]*b.d[66] + a.d[43]*b.d[67] + a.d[52]*b.d[68] + a.d[61]*b.d[69] + a.d[70]*b.d[70] + a.d[79]*b.d[71];
    c.d[71] = a.d[7]*b.d[72] + a.d[16]*b.d[73] + a.d[25]*b.d[74] + a.d[34]*b.d[75] + a.d[43]*b.d[76] + a.d[52]*b.d[77] + a.d[61]*b.d[78] + a.d[70]*b.d[79] + a.d[79]*b.d[80];

    c.d[72] = a.d[8]*b.d[0] + a.d[17]*b.d[1] + a.d[26]*b.d[2] + a.d[35]*b.d[3] + a.d[44]*b.d[4] + a.d[53]*b.d[5] + a.d[62]*b.d[6] + a.d[71]*b.d[7] + a.d[80]*b.d[8];
    c.d[73] = a.d[8]*b.d[9] + a.d[17]*b.d[10] + a.d[26]*b.d[11] + a.d[35]*b.d[12] + a.d[44]*b.d[13] + a.d[53]*b.d[14] + a.d[62]*b.d[15] + a.d[71]*b.d[16] + a.d[80]*b.d[17];
    c.d[74] = a.d[8]*b.d[18] + a.d[17]*b.d[19] + a.d[26]*b.d[20] + a.d[35]*b.d[21] + a.d[44]*b.d[22] + a.d[53]*b.d[23] + a.d[62]*b.d[24] + a.d[71]*b.d[25] + a.d[80]*b.d[26];
    c.d[75] = a.d[8]*b.d[27] + a.d[17]*b.d[28] + a.d[26]*b.d[29] + a.d[35]*b.d[30] + a.d[44]*b.d[31] + a.d[53]*b.d[32] + a.d[62]*b.d[33] + a.d[71]*b.d[34] + a.d[80]*b.d[35];
    c.d[76] = a.d[8]*b.d[36] + a.d[17]*b.d[37] + a.d[26]*b.d[38] + a.d[35]*b.d[39] + a.d[44]*b.d[40] + a.d[53]*b.d[41] + a.d[62]*b.d[42] + a.d[71]*b.d[43] + a.d[80]*b.d[44];
    c.d[77] = a.d[8]*b.d[45] + a.d[17]*b.d[46] + a.d[26]*b.d[47] + a.d[35]*b.d[48] + a.d[44]*b.d[49] + a.d[53]*b.d[50] + a.d[62]*b.d[51] + a.d[71]*b.d[52] + a.d[80]*b.d[53];
    c.d[78] = a.d[8]*b.d[54] + a.d[17]*b.d[55] + a.d[26]*b.d[56] + a.d[35]*b.d[57] + a.d[44]*b.d[58] + a.d[53]*b.d[59] + a.d[62]*b.d[60] + a.d[71]*b.d[61] + a.d[80]*b.d[62];
    c.d[79] = a.d[8]*b.d[63] + a.d[17]*b.d[64] + a.d[26]*b.d[65] + a.d[35]*b.d[66] + a.d[44]*b.d[67] + a.d[53]*b.d[68] + a.d[62]*b.d[69] + a.d[71]*b.d[70] + a.d[80]*b.d[71];
    c.d[80] = a.d[8]*b.d[72] + a.d[17]*b.d[73] + a.d[26]*b.d[74] + a.d[35]*b.d[75] + a.d[44]*b.d[76] + a.d[53]*b.d[77] + a.d[62]*b.d[78] + a.d[71]*b.d[79] + a.d[80]*b.d[80];

    return c;
    
}

//-----------------------------------------------------------------------------
// (a ddot b)_ijkl = a_ijmn b_mnkl where b is super-symmetric
inline tens4d ddot(const tens4d& a, const tens4ds& b)
{
    tens4d c;
    
    c.d[0] = a.d[0]*b.d[0] + a.d[9]*b.d[1] + a.d[18]*b.d[3] + a.d[27]*b.d[6] + a.d[54]*b.d[6] + a.d[36]*b.d[10] + a.d[63]*b.d[10] + a.d[45]*b.d[15] + a.d[72]*b.d[15];
    c.d[1] = a.d[0]*b.d[1] + a.d[9]*b.d[2] + a.d[18]*b.d[4] + a.d[27]*b.d[7] + a.d[54]*b.d[7] + a.d[36]*b.d[11] + a.d[63]*b.d[11] + a.d[45]*b.d[16] + a.d[72]*b.d[16];
    c.d[2] = a.d[0]*b.d[3] + a.d[9]*b.d[4] + a.d[18]*b.d[5] + a.d[27]*b.d[8] + a.d[54]*b.d[8] + a.d[36]*b.d[12] + a.d[63]*b.d[12] + a.d[45]*b.d[17] + a.d[72]*b.d[17];
    c.d[3] = a.d[0]*b.d[6] + a.d[9]*b.d[7] + a.d[18]*b.d[8] + a.d[27]*b.d[9] + a.d[54]*b.d[9] + a.d[36]*b.d[13] + a.d[63]*b.d[13] + a.d[45]*b.d[18] + a.d[72]*b.d[18];
    c.d[4] = a.d[0]*b.d[10] + a.d[9]*b.d[11] + a.d[18]*b.d[12] + a.d[27]*b.d[13] + a.d[54]*b.d[13] + a.d[36]*b.d[14] + a.d[63]*b.d[14] + a.d[45]*b.d[19] + a.d[72]*b.d[19];
    c.d[5] = a.d[0]*b.d[15] + a.d[9]*b.d[16] + a.d[18]*b.d[17] + a.d[27]*b.d[18] + a.d[54]*b.d[18] + a.d[36]*b.d[19] + a.d[63]*b.d[19] + a.d[45]*b.d[20] + a.d[72]*b.d[20];
    c.d[6] = a.d[0]*b.d[6] + a.d[9]*b.d[7] + a.d[18]*b.d[8] + a.d[27]*b.d[9] + a.d[54]*b.d[9] + a.d[36]*b.d[13] + a.d[63]*b.d[13] + a.d[45]*b.d[18] + a.d[72]*b.d[18];
    c.d[7] = a.d[0]*b.d[10] + a.d[9]*b.d[11] + a.d[18]*b.d[12] + a.d[27]*b.d[13] + a.d[54]*b.d[13] + a.d[36]*b.d[14] + a.d[63]*b.d[14] + a.d[45]*b.d[19] + a.d[72]*b.d[19];
    c.d[8] = a.d[0]*b.d[15] + a.d[9]*b.d[16] + a.d[18]*b.d[17] + a.d[27]*b.d[18] + a.d[54]*b.d[18] + a.d[36]*b.d[19] + a.d[63]*b.d[19] + a.d[45]*b.d[20] + a.d[72]*b.d[20];

    c.d[9] = a.d[1]*b.d[0] + a.d[10]*b.d[1] + a.d[19]*b.d[3] + a.d[28]*b.d[6] + a.d[55]*b.d[6] + a.d[37]*b.d[10] + a.d[64]*b.d[10] + a.d[46]*b.d[15] + a.d[73]*b.d[15];
    c.d[10] = a.d[1]*b.d[1] + a.d[10]*b.d[2] + a.d[19]*b.d[4] + a.d[28]*b.d[7] + a.d[55]*b.d[7] + a.d[37]*b.d[11] + a.d[64]*b.d[11] + a.d[46]*b.d[16] + a.d[73]*b.d[16];
    c.d[11] = a.d[1]*b.d[3] + a.d[10]*b.d[4] + a.d[19]*b.d[5] + a.d[28]*b.d[8] + a.d[55]*b.d[8] + a.d[37]*b.d[12] + a.d[64]*b.d[12] + a.d[46]*b.d[17] + a.d[73]*b.d[17];
    c.d[12] = a.d[1]*b.d[6] + a.d[10]*b.d[7] + a.d[19]*b.d[8] + a.d[28]*b.d[9] + a.d[55]*b.d[9] + a.d[37]*b.d[13] + a.d[64]*b.d[13] + a.d[46]*b.d[18] + a.d[73]*b.d[18];
    c.d[13] = a.d[1]*b.d[10] + a.d[10]*b.d[11] + a.d[19]*b.d[12] + a.d[28]*b.d[13] + a.d[55]*b.d[13] + a.d[37]*b.d[14] + a.d[64]*b.d[14] + a.d[46]*b.d[19] + a.d[73]*b.d[19];
    c.d[14] = a.d[1]*b.d[15] + a.d[10]*b.d[16] + a.d[19]*b.d[17] + a.d[28]*b.d[18] + a.d[55]*b.d[18] + a.d[37]*b.d[19] + a.d[64]*b.d[19] + a.d[46]*b.d[20] + a.d[73]*b.d[20];
    c.d[15] = a.d[1]*b.d[6] + a.d[10]*b.d[7] + a.d[19]*b.d[8] + a.d[28]*b.d[9] + a.d[55]*b.d[9] + a.d[37]*b.d[13] + a.d[64]*b.d[13] + a.d[46]*b.d[18] + a.d[73]*b.d[18];
    c.d[16] = a.d[1]*b.d[10] + a.d[10]*b.d[11] + a.d[19]*b.d[12] + a.d[28]*b.d[13] + a.d[55]*b.d[13] + a.d[37]*b.d[14] + a.d[64]*b.d[14] + a.d[46]*b.d[19] + a.d[73]*b.d[19];
    c.d[17] = a.d[1]*b.d[15] + a.d[10]*b.d[16] + a.d[19]*b.d[17] + a.d[28]*b.d[18] + a.d[55]*b.d[18] + a.d[37]*b.d[19] + a.d[64]*b.d[19] + a.d[46]*b.d[20] + a.d[73]*b.d[20];

    c.d[18] = a.d[2]*b.d[0] + a.d[11]*b.d[1] + a.d[20]*b.d[3] + a.d[29]*b.d[6] + a.d[56]*b.d[6] + a.d[38]*b.d[10] + a.d[65]*b.d[10] + a.d[47]*b.d[15] + a.d[74]*b.d[15];
    c.d[19] = a.d[2]*b.d[1] + a.d[11]*b.d[2] + a.d[20]*b.d[4] + a.d[29]*b.d[7] + a.d[56]*b.d[7] + a.d[38]*b.d[11] + a.d[65]*b.d[11] + a.d[47]*b.d[16] + a.d[74]*b.d[16];
    c.d[20] = a.d[2]*b.d[3] + a.d[11]*b.d[4] + a.d[20]*b.d[5] + a.d[29]*b.d[8] + a.d[56]*b.d[8] + a.d[38]*b.d[12] + a.d[65]*b.d[12] + a.d[47]*b.d[17] + a.d[74]*b.d[17];
    c.d[21] = a.d[2]*b.d[6] + a.d[11]*b.d[7] + a.d[20]*b.d[8] + a.d[29]*b.d[9] + a.d[56]*b.d[9] + a.d[38]*b.d[13] + a.d[65]*b.d[13] + a.d[47]*b.d[18] + a.d[74]*b.d[18];
    c.d[22] = a.d[2]*b.d[10] + a.d[11]*b.d[11] + a.d[20]*b.d[12] + a.d[29]*b.d[13] + a.d[56]*b.d[13] + a.d[38]*b.d[14] + a.d[65]*b.d[14] + a.d[47]*b.d[19] + a.d[74]*b.d[19];
    c.d[23] = a.d[2]*b.d[15] + a.d[11]*b.d[16] + a.d[20]*b.d[17] + a.d[29]*b.d[18] + a.d[56]*b.d[18] + a.d[38]*b.d[19] + a.d[65]*b.d[19] + a.d[47]*b.d[20] + a.d[74]*b.d[20];
    c.d[24] = a.d[2]*b.d[6] + a.d[11]*b.d[7] + a.d[20]*b.d[8] + a.d[29]*b.d[9] + a.d[56]*b.d[9] + a.d[38]*b.d[13] + a.d[65]*b.d[13] + a.d[47]*b.d[18] + a.d[74]*b.d[18];
    c.d[25] = a.d[2]*b.d[10] + a.d[11]*b.d[11] + a.d[20]*b.d[12] + a.d[29]*b.d[13] + a.d[56]*b.d[13] + a.d[38]*b.d[14] + a.d[65]*b.d[14] + a.d[47]*b.d[19] + a.d[74]*b.d[19];
    c.d[26] = a.d[2]*b.d[15] + a.d[11]*b.d[16] + a.d[20]*b.d[17] + a.d[29]*b.d[18] + a.d[56]*b.d[18] + a.d[38]*b.d[19] + a.d[65]*b.d[19] + a.d[47]*b.d[20] + a.d[74]*b.d[20];

    c.d[27] = a.d[3]*b.d[0] + a.d[12]*b.d[1] + a.d[21]*b.d[3] + a.d[30]*b.d[6] + a.d[57]*b.d[6] + a.d[39]*b.d[10] + a.d[66]*b.d[10] + a.d[48]*b.d[15] + a.d[75]*b.d[15];
    c.d[28] = a.d[3]*b.d[1] + a.d[12]*b.d[2] + a.d[21]*b.d[4] + a.d[30]*b.d[7] + a.d[57]*b.d[7] + a.d[39]*b.d[11] + a.d[66]*b.d[11] + a.d[48]*b.d[16] + a.d[75]*b.d[16];
    c.d[29] = a.d[3]*b.d[3] + a.d[12]*b.d[4] + a.d[21]*b.d[5] + a.d[30]*b.d[8] + a.d[57]*b.d[8] + a.d[39]*b.d[12] + a.d[66]*b.d[12] + a.d[48]*b.d[17] + a.d[75]*b.d[17];
    c.d[30] = a.d[3]*b.d[6] + a.d[12]*b.d[7] + a.d[21]*b.d[8] + a.d[30]*b.d[9] + a.d[57]*b.d[9] + a.d[39]*b.d[13] + a.d[66]*b.d[13] + a.d[48]*b.d[18] + a.d[75]*b.d[18];
    c.d[31] = a.d[3]*b.d[10] + a.d[12]*b.d[11] + a.d[21]*b.d[12] + a.d[30]*b.d[13] + a.d[57]*b.d[13] + a.d[39]*b.d[14] + a.d[66]*b.d[14] + a.d[48]*b.d[19] + a.d[75]*b.d[19];
    c.d[32] = a.d[3]*b.d[15] + a.d[12]*b.d[16] + a.d[21]*b.d[17] + a.d[30]*b.d[18] + a.d[57]*b.d[18] + a.d[39]*b.d[19] + a.d[66]*b.d[19] + a.d[48]*b.d[20] + a.d[75]*b.d[20];
    c.d[33] = a.d[3]*b.d[6] + a.d[12]*b.d[7] + a.d[21]*b.d[8] + a.d[30]*b.d[9] + a.d[57]*b.d[9] + a.d[39]*b.d[13] + a.d[66]*b.d[13] + a.d[48]*b.d[18] + a.d[75]*b.d[18];
    c.d[34] = a.d[3]*b.d[10] + a.d[12]*b.d[11] + a.d[21]*b.d[12] + a.d[30]*b.d[13] + a.d[57]*b.d[13] + a.d[39]*b.d[14] + a.d[66]*b.d[14] + a.d[48]*b.d[19] + a.d[75]*b.d[19];
    c.d[35] = a.d[3]*b.d[15] + a.d[12]*b.d[16] + a.d[21]*b.d[17] + a.d[30]*b.d[18] + a.d[57]*b.d[18] + a.d[39]*b.d[19] + a.d[66]*b.d[19] + a.d[48]*b.d[20] + a.d[75]*b.d[20];

    c.d[36] = a.d[4]*b.d[0] + a.d[13]*b.d[1] + a.d[22]*b.d[3] + a.d[31]*b.d[6] + a.d[58]*b.d[6] + a.d[40]*b.d[10] + a.d[67]*b.d[10] + a.d[49]*b.d[15] + a.d[76]*b.d[15];
    c.d[37] = a.d[4]*b.d[1] + a.d[13]*b.d[2] + a.d[22]*b.d[4] + a.d[31]*b.d[7] + a.d[58]*b.d[7] + a.d[40]*b.d[11] + a.d[67]*b.d[11] + a.d[49]*b.d[16] + a.d[76]*b.d[16];
    c.d[38] = a.d[4]*b.d[3] + a.d[13]*b.d[4] + a.d[22]*b.d[5] + a.d[31]*b.d[8] + a.d[58]*b.d[8] + a.d[40]*b.d[12] + a.d[67]*b.d[12] + a.d[49]*b.d[17] + a.d[76]*b.d[17];
    c.d[39] = a.d[4]*b.d[6] + a.d[13]*b.d[7] + a.d[22]*b.d[8] + a.d[31]*b.d[9] + a.d[58]*b.d[9] + a.d[40]*b.d[13] + a.d[67]*b.d[13] + a.d[49]*b.d[18] + a.d[76]*b.d[18];
    c.d[40] = a.d[4]*b.d[10] + a.d[13]*b.d[11] + a.d[22]*b.d[12] + a.d[31]*b.d[13] + a.d[58]*b.d[13] + a.d[40]*b.d[14] + a.d[67]*b.d[14] + a.d[49]*b.d[19] + a.d[76]*b.d[19];
    c.d[41] = a.d[4]*b.d[15] + a.d[13]*b.d[16] + a.d[22]*b.d[17] + a.d[31]*b.d[18] + a.d[58]*b.d[18] + a.d[40]*b.d[19] + a.d[67]*b.d[19] + a.d[49]*b.d[20] + a.d[76]*b.d[20];
    c.d[42] = a.d[4]*b.d[6] + a.d[13]*b.d[7] + a.d[22]*b.d[8] + a.d[31]*b.d[9] + a.d[58]*b.d[9] + a.d[40]*b.d[13] + a.d[67]*b.d[13] + a.d[49]*b.d[18] + a.d[76]*b.d[18];
    c.d[43] = a.d[4]*b.d[10] + a.d[13]*b.d[11] + a.d[22]*b.d[12] + a.d[31]*b.d[13] + a.d[58]*b.d[13] + a.d[40]*b.d[14] + a.d[67]*b.d[14] + a.d[49]*b.d[19] + a.d[76]*b.d[19];
    c.d[44] = a.d[4]*b.d[15] + a.d[13]*b.d[16] + a.d[22]*b.d[17] + a.d[31]*b.d[18] + a.d[58]*b.d[18] + a.d[40]*b.d[19] + a.d[67]*b.d[19] + a.d[49]*b.d[20] + a.d[76]*b.d[20];

    c.d[45] = a.d[5]*b.d[0] + a.d[14]*b.d[1] + a.d[23]*b.d[3] + a.d[32]*b.d[6] + a.d[59]*b.d[6] + a.d[41]*b.d[10] + a.d[68]*b.d[10] + a.d[50]*b.d[15] + a.d[77]*b.d[15];
    c.d[46] = a.d[5]*b.d[1] + a.d[14]*b.d[2] + a.d[23]*b.d[4] + a.d[32]*b.d[7] + a.d[59]*b.d[7] + a.d[41]*b.d[11] + a.d[68]*b.d[11] + a.d[50]*b.d[16] + a.d[77]*b.d[16];
    c.d[47] = a.d[5]*b.d[3] + a.d[14]*b.d[4] + a.d[23]*b.d[5] + a.d[32]*b.d[8] + a.d[59]*b.d[8] + a.d[41]*b.d[12] + a.d[68]*b.d[12] + a.d[50]*b.d[17] + a.d[77]*b.d[17];
    c.d[48] = a.d[5]*b.d[6] + a.d[14]*b.d[7] + a.d[23]*b.d[8] + a.d[32]*b.d[9] + a.d[59]*b.d[9] + a.d[41]*b.d[13] + a.d[68]*b.d[13] + a.d[50]*b.d[18] + a.d[77]*b.d[18];
    c.d[49] = a.d[5]*b.d[10] + a.d[14]*b.d[11] + a.d[23]*b.d[12] + a.d[32]*b.d[13] + a.d[59]*b.d[13] + a.d[41]*b.d[14] + a.d[68]*b.d[14] + a.d[50]*b.d[19] + a.d[77]*b.d[19];
    c.d[50] = a.d[5]*b.d[15] + a.d[14]*b.d[16] + a.d[23]*b.d[17] + a.d[32]*b.d[18] + a.d[59]*b.d[18] + a.d[41]*b.d[19] + a.d[68]*b.d[19] + a.d[50]*b.d[20] + a.d[77]*b.d[20];
    c.d[51] = a.d[5]*b.d[6] + a.d[14]*b.d[7] + a.d[23]*b.d[8] + a.d[32]*b.d[9] + a.d[59]*b.d[9] + a.d[41]*b.d[13] + a.d[68]*b.d[13] + a.d[50]*b.d[18] + a.d[77]*b.d[18];
    c.d[52] = a.d[5]*b.d[10] + a.d[14]*b.d[11] + a.d[23]*b.d[12] + a.d[32]*b.d[13] + a.d[59]*b.d[13] + a.d[41]*b.d[14] + a.d[68]*b.d[14] + a.d[50]*b.d[19] + a.d[77]*b.d[19];
    c.d[53] = a.d[5]*b.d[15] + a.d[14]*b.d[16] + a.d[23]*b.d[17] + a.d[32]*b.d[18] + a.d[59]*b.d[18] + a.d[41]*b.d[19] + a.d[68]*b.d[19] + a.d[50]*b.d[20] + a.d[77]*b.d[20];

    c.d[54] = a.d[6]*b.d[0] + a.d[15]*b.d[1] + a.d[24]*b.d[3] + a.d[33]*b.d[6] + a.d[60]*b.d[6] + a.d[42]*b.d[10] + a.d[69]*b.d[10] + a.d[51]*b.d[15] + a.d[78]*b.d[15];
    c.d[55] = a.d[6]*b.d[1] + a.d[15]*b.d[2] + a.d[24]*b.d[4] + a.d[33]*b.d[7] + a.d[60]*b.d[7] + a.d[42]*b.d[11] + a.d[69]*b.d[11] + a.d[51]*b.d[16] + a.d[78]*b.d[16];
    c.d[56] = a.d[6]*b.d[3] + a.d[15]*b.d[4] + a.d[24]*b.d[5] + a.d[33]*b.d[8] + a.d[60]*b.d[8] + a.d[42]*b.d[12] + a.d[69]*b.d[12] + a.d[51]*b.d[17] + a.d[78]*b.d[17];
    c.d[57] = a.d[6]*b.d[6] + a.d[15]*b.d[7] + a.d[24]*b.d[8] + a.d[33]*b.d[9] + a.d[60]*b.d[9] + a.d[42]*b.d[13] + a.d[69]*b.d[13] + a.d[51]*b.d[18] + a.d[78]*b.d[18];
    c.d[58] = a.d[6]*b.d[10] + a.d[15]*b.d[11] + a.d[24]*b.d[12] + a.d[33]*b.d[13] + a.d[60]*b.d[13] + a.d[42]*b.d[14] + a.d[69]*b.d[14] + a.d[51]*b.d[19] + a.d[78]*b.d[19];
    c.d[59] = a.d[6]*b.d[15] + a.d[15]*b.d[16] + a.d[24]*b.d[17] + a.d[33]*b.d[18] + a.d[60]*b.d[18] + a.d[42]*b.d[19] + a.d[69]*b.d[19] + a.d[51]*b.d[20] + a.d[78]*b.d[20];
    c.d[60] = a.d[6]*b.d[6] + a.d[15]*b.d[7] + a.d[24]*b.d[8] + a.d[33]*b.d[9] + a.d[60]*b.d[9] + a.d[42]*b.d[13] + a.d[69]*b.d[13] + a.d[51]*b.d[18] + a.d[78]*b.d[18];
    c.d[61] = a.d[6]*b.d[10] + a.d[15]*b.d[11] + a.d[24]*b.d[12] + a.d[33]*b.d[13] + a.d[60]*b.d[13] + a.d[42]*b.d[14] + a.d[69]*b.d[14] + a.d[51]*b.d[19] + a.d[78]*b.d[19];
    c.d[62] = a.d[6]*b.d[15] + a.d[15]*b.d[16] + a.d[24]*b.d[17] + a.d[33]*b.d[18] + a.d[60]*b.d[18] + a.d[42]*b.d[19] + a.d[69]*b.d[19] + a.d[51]*b.d[20] + a.d[78]*b.d[20];

    c.d[63] = a.d[7]*b.d[0] + a.d[16]*b.d[1] + a.d[25]*b.d[3] + a.d[34]*b.d[6] + a.d[61]*b.d[6] + a.d[43]*b.d[10] + a.d[70]*b.d[10] + a.d[52]*b.d[15] + a.d[79]*b.d[15];
    c.d[64] = a.d[7]*b.d[1] + a.d[16]*b.d[2] + a.d[25]*b.d[4] + a.d[34]*b.d[7] + a.d[61]*b.d[7] + a.d[43]*b.d[11] + a.d[70]*b.d[11] + a.d[52]*b.d[16] + a.d[79]*b.d[16];
    c.d[65] = a.d[7]*b.d[3] + a.d[16]*b.d[4] + a.d[25]*b.d[5] + a.d[34]*b.d[8] + a.d[61]*b.d[8] + a.d[43]*b.d[12] + a.d[70]*b.d[12] + a.d[52]*b.d[17] + a.d[79]*b.d[17];
    c.d[66] = a.d[7]*b.d[6] + a.d[16]*b.d[7] + a.d[25]*b.d[8] + a.d[34]*b.d[9] + a.d[61]*b.d[9] + a.d[43]*b.d[13] + a.d[70]*b.d[13] + a.d[52]*b.d[18] + a.d[79]*b.d[18];
    c.d[67] = a.d[7]*b.d[10] + a.d[16]*b.d[11] + a.d[25]*b.d[12] + a.d[34]*b.d[13] + a.d[61]*b.d[13] + a.d[43]*b.d[14] + a.d[70]*b.d[14] + a.d[52]*b.d[19] + a.d[79]*b.d[19];
    c.d[68] = a.d[7]*b.d[15] + a.d[16]*b.d[16] + a.d[25]*b.d[17] + a.d[34]*b.d[18] + a.d[61]*b.d[18] + a.d[43]*b.d[19] + a.d[70]*b.d[19] + a.d[52]*b.d[20] + a.d[79]*b.d[20];
    c.d[69] = a.d[7]*b.d[6] + a.d[16]*b.d[7] + a.d[25]*b.d[8] + a.d[34]*b.d[9] + a.d[61]*b.d[9] + a.d[43]*b.d[13] + a.d[70]*b.d[13] + a.d[52]*b.d[18] + a.d[79]*b.d[18];
    c.d[70] = a.d[7]*b.d[10] + a.d[16]*b.d[11] + a.d[25]*b.d[12] + a.d[34]*b.d[13] + a.d[61]*b.d[13] + a.d[43]*b.d[14] + a.d[70]*b.d[14] + a.d[52]*b.d[19] + a.d[79]*b.d[19];
    c.d[71] = a.d[7]*b.d[15] + a.d[16]*b.d[16] + a.d[25]*b.d[17] + a.d[34]*b.d[18] + a.d[61]*b.d[18] + a.d[43]*b.d[19] + a.d[70]*b.d[19] + a.d[52]*b.d[20] + a.d[79]*b.d[20];

    c.d[72] = a.d[8]*b.d[0] + a.d[17]*b.d[1] + a.d[26]*b.d[3] + a.d[35]*b.d[6] + a.d[62]*b.d[6] + a.d[44]*b.d[10] + a.d[71]*b.d[10] + a.d[53]*b.d[15] + a.d[80]*b.d[15];
    c.d[73] = a.d[8]*b.d[1] + a.d[17]*b.d[2] + a.d[26]*b.d[4] + a.d[35]*b.d[7] + a.d[62]*b.d[7] + a.d[44]*b.d[11] + a.d[71]*b.d[11] + a.d[53]*b.d[16] + a.d[80]*b.d[16];
    c.d[74] = a.d[8]*b.d[3] + a.d[17]*b.d[4] + a.d[26]*b.d[5] + a.d[35]*b.d[8] + a.d[62]*b.d[8] + a.d[44]*b.d[12] + a.d[71]*b.d[12] + a.d[53]*b.d[17] + a.d[80]*b.d[17];
    c.d[75] = a.d[8]*b.d[6] + a.d[17]*b.d[7] + a.d[26]*b.d[8] + a.d[35]*b.d[9] + a.d[62]*b.d[9] + a.d[44]*b.d[13] + a.d[71]*b.d[13] + a.d[53]*b.d[18] + a.d[80]*b.d[18];
    c.d[76] = a.d[8]*b.d[10] + a.d[17]*b.d[11] + a.d[26]*b.d[12] + a.d[35]*b.d[13] + a.d[62]*b.d[13] + a.d[44]*b.d[14] + a.d[71]*b.d[14] + a.d[53]*b.d[19] + a.d[80]*b.d[19];
    c.d[77] = a.d[8]*b.d[15] + a.d[17]*b.d[16] + a.d[26]*b.d[17] + a.d[35]*b.d[18] + a.d[62]*b.d[18] + a.d[44]*b.d[19] + a.d[71]*b.d[19] + a.d[53]*b.d[20] + a.d[80]*b.d[20];
    c.d[78] = a.d[8]*b.d[6] + a.d[17]*b.d[7] + a.d[26]*b.d[8] + a.d[35]*b.d[9] + a.d[62]*b.d[9] + a.d[44]*b.d[13] + a.d[71]*b.d[13] + a.d[53]*b.d[18] + a.d[80]*b.d[18];
    c.d[79] = a.d[8]*b.d[10] + a.d[17]*b.d[11] + a.d[26]*b.d[12] + a.d[35]*b.d[13] + a.d[62]*b.d[13] + a.d[44]*b.d[14] + a.d[71]*b.d[14] + a.d[53]*b.d[19] + a.d[80]*b.d[19];
    c.d[80] = a.d[8]*b.d[15] + a.d[17]*b.d[16] + a.d[26]*b.d[17] + a.d[35]*b.d[18] + a.d[62]*b.d[18] + a.d[44]*b.d[19] + a.d[71]*b.d[19] + a.d[53]*b.d[20] + a.d[80]*b.d[20];

    return c;
    
}

//-----------------------------------------------------------------------------
// vdotTdotv_jk = a_i T_ijkl b_l
inline mat3d vdotTdotv(const vec3d& a, const tens4d& T, const vec3d& b)
{
    return mat3d(a.x*b.x*T.d[0] + a.z*b.x*T.d[5] + a.y*b.x*T.d[6] + a.x*b.y*T.d[27] + a.z*b.y*T.d[32] + a.y*b.y*T.d[33] + a.x*b.z*T.d[72] + a.z*b.z*T.d[77] + a.y*b.z*T.d[78],
    a.y*b.x*T.d[1] + a.x*b.x*T.d[3] + a.z*b.x*T.d[7] + a.y*b.y*T.d[28] + a.x*b.y*T.d[30] + a.z*b.y*T.d[34] + a.y*b.z*T.d[73] + a.x*b.z*T.d[75] + a.z*b.z*T.d[79],
    a.z*b.x*T.d[2] + a.y*b.x*T.d[4] + a.x*b.x*T.d[8] + a.z*b.y*T.d[29] + a.y*b.y*T.d[31] + a.x*b.y*T.d[35] + a.z*b.z*T.d[74] + a.y*b.z*T.d[76] + a.x*b.z*T.d[80],
    a.x*b.y*T.d[9] + a.z*b.y*T.d[14] + a.y*b.y*T.d[15] + a.x*b.z*T.d[36] + a.z*b.z*T.d[41] + a.y*b.z*T.d[42] + a.x*b.x*T.d[54] + a.z*b.x*T.d[59] + a.y*b.x*T.d[60],
    a.y*b.y*T.d[10] + a.x*b.y*T.d[12] + a.z*b.y*T.d[16] + a.y*b.z*T.d[37] + a.x*b.z*T.d[39] + a.z*b.z*T.d[43] + a.y*b.x*T.d[55] + a.x*b.x*T.d[57] + a.z*b.x*T.d[61],
    a.z*b.y*T.d[11] + a.y*b.y*T.d[13] + a.x*b.y*T.d[17] + a.z*b.z*T.d[38] + a.y*b.z*T.d[40] + a.x*b.z*T.d[44] + a.z*b.x*T.d[56] + a.y*b.x*T.d[58] + a.x*b.x*T.d[62],
    a.x*b.z*T.d[18] + a.z*b.z*T.d[23] + a.y*b.z*T.d[24] + a.x*b.x*T.d[45] + a.z*b.x*T.d[50] + a.y*b.x*T.d[51] + a.x*b.y*T.d[63] + a.z*b.y*T.d[68] + a.y*b.y*T.d[69],
    a.y*b.z*T.d[19] + a.x*b.z*T.d[21] + a.z*b.z*T.d[25] + a.y*b.x*T.d[46] + a.x*b.x*T.d[48] + a.z*b.x*T.d[52] + a.y*b.y*T.d[64] + a.x*b.y*T.d[66] + a.z*b.y*T.d[70],
    a.z*b.z*T.d[20] + a.y*b.z*T.d[22] + a.x*b.z*T.d[26] + a.z*b.x*T.d[47] + a.y*b.x*T.d[49] + a.x*b.x*T.d[53] + a.z*b.y*T.d[65] + a.y*b.y*T.d[67] + a.x*b.y*T.d[71]);
}

//-----------------------------------------------------------------------------
// inverse
inline tens4d tens4d::inverse() const
{
    matrix c(9,9);
    
    // populate c
    c(0,0) = d[ 0]; c(0,1) = d[ 9]; c(0,2) = d[18]; c(0,3) = d[27]; c(0,4) = d[36]; c(0,5) = d[45]; c(0,6) = d[54]; c(0,7) = d[63]; c(0,8) = d[72];
    c(1,0) = d[ 1]; c(1,1) = d[10]; c(1,2) = d[19]; c(1,3) = d[28]; c(1,4) = d[37]; c(1,5) = d[46]; c(1,6) = d[55]; c(1,7) = d[64]; c(1,8) = d[73];
    c(2,0) = d[ 2]; c(2,1) = d[11]; c(2,2) = d[20]; c(2,3) = d[29]; c(2,4) = d[38]; c(2,5) = d[47]; c(2,6) = d[56]; c(2,7) = d[65]; c(2,8) = d[74];
    c(3,0) = d[ 3]; c(3,1) = d[12]; c(3,2) = d[21]; c(3,3) = d[30]; c(3,4) = d[39]; c(3,5) = d[48]; c(3,6) = d[57]; c(3,7) = d[66]; c(3,8) = d[75];
    c(4,0) = d[ 4]; c(4,1) = d[13]; c(4,2) = d[22]; c(4,3) = d[31]; c(4,4) = d[40]; c(4,5) = d[49]; c(4,6) = d[58]; c(4,7) = d[67]; c(4,8) = d[76];
    c(5,0) = d[ 5]; c(5,1) = d[14]; c(5,2) = d[23]; c(5,3) = d[32]; c(5,4) = d[41]; c(5,5) = d[50]; c(5,6) = d[59]; c(5,7) = d[68]; c(5,8) = d[77];
    c(6,0) = d[ 6]; c(6,1) = d[15]; c(6,2) = d[24]; c(6,3) = d[33]; c(6,4) = d[42]; c(6,5) = d[51]; c(6,6) = d[60]; c(6,7) = d[69]; c(6,8) = d[78];
    c(7,0) = d[ 7]; c(7,1) = d[16]; c(7,2) = d[25]; c(7,3) = d[34]; c(7,4) = d[43]; c(7,5) = d[52]; c(7,6) = d[61]; c(7,7) = d[70]; c(7,8) = d[79];
    c(8,0) = d[ 8]; c(8,1) = d[17]; c(8,2) = d[26]; c(8,3) = d[35]; c(8,4) = d[44]; c(8,5) = d[53]; c(8,6) = d[62]; c(8,7) = d[71]; c(8,8) = d[80];

    // invert c
    matrix s = c.inverse();
    
    // return inverse
    tens4d S;
    S.d[ 0] = s(0,0); S.d[ 9] = s(0,1); S.d[18] = s(0,2); S.d[27] = s(0,3); S.d[36] = s(0,4); S.d[45] = s(0,5); S.d[54] = s(0,6); S.d[63] = s(0,7); S.d[72] = s(0,8);
    S.d[ 1] = s(1,0); S.d[10] = s(1,1); S.d[19] = s(1,2); S.d[28] = s(1,3); S.d[37] = s(1,4); S.d[46] = s(1,5); S.d[55] = s(1,6); S.d[64] = s(1,7); S.d[73] = s(1,8);
    S.d[ 2] = s(2,0); S.d[11] = s(2,1); S.d[20] = s(2,2); S.d[29] = s(2,3); S.d[38] = s(2,4); S.d[47] = s(2,5); S.d[56] = s(2,6); S.d[65] = s(2,7); S.d[74] = s(2,8);
    S.d[ 3] = s(3,0); S.d[12] = s(3,1); S.d[21] = s(3,2); S.d[30] = s(3,3); S.d[39] = s(3,4); S.d[48] = s(3,5); S.d[57] = s(3,6); S.d[66] = s(3,7); S.d[75] = s(3,8);
    S.d[ 4] = s(4,0); S.d[13] = s(4,1); S.d[22] = s(4,2); S.d[31] = s(4,3); S.d[40] = s(4,4); S.d[49] = s(4,5); S.d[58] = s(4,6); S.d[67] = s(4,7); S.d[76] = s(4,8);
    S.d[ 5] = s(5,0); S.d[14] = s(5,1); S.d[23] = s(5,2); S.d[32] = s(5,3); S.d[41] = s(5,4); S.d[50] = s(5,5); S.d[59] = s(5,6); S.d[68] = s(5,7); S.d[77] = s(5,8);
    S.d[ 6] = s(6,0); S.d[15] = s(6,1); S.d[24] = s(6,2); S.d[33] = s(6,3); S.d[42] = s(6,4); S.d[51] = s(6,5); S.d[60] = s(6,6); S.d[69] = s(6,7); S.d[78] = s(6,8);
    S.d[ 7] = s(7,0); S.d[16] = s(7,1); S.d[25] = s(7,2); S.d[34] = s(7,3); S.d[43] = s(7,4); S.d[52] = s(7,5); S.d[61] = s(7,6); S.d[70] = s(7,7); S.d[79] = s(7,8);
    S.d[ 8] = s(8,0); S.d[17] = s(8,1); S.d[26] = s(8,2); S.d[35] = s(8,3); S.d[44] = s(8,4); S.d[53] = s(8,5); S.d[62] = s(8,6); S.d[71] = s(8,7); S.d[80] = s(8,8);

    return S;
}

